# Does small-scale turbulence matter for ice growth in mixed-phase clouds?

G. Sarnitsky, G. Sardina, G. Svensson, Alain Pumir, F. Hoffmann, B. Mehlig

► To cite this version:

G. Sarnitsky, G. Sardina, G. Svensson, Alain Pumir, F. Hoffmann, et al. Does small-scale turbulence matter for ice growth in mixed-phase clouds?. Physical Review Fluids, 2025, 10 (5), pp.053803. 10.1103/PhysRevFluids.10.053803. hal-05375082

HAL Id: hal-05375082

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

Submitted on 20 Nov 2025

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 Attribution 4.0 International License
Does small-scale turbulence matter for ice growth in mixed-phase clouds?

G. Sarnitsky,¹ G. Sardina,² G. Svensson,³,⁴ A. Pumir,⁵ F. Hoffmann,⁶ and B. Mehlig¹, *

¹Department of Physics, Gothenburg University, SE-41296 Gothenburg, Sweden

²Department of Mechanics and Maritime Sciences,
Chalmers University of Technology, SE-41296 Gothenburg, Sweden

³Department of Meteorology and Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden

⁴Department of Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden

⁵Univ. Lyon, ENS de Lyon, Univ. Claude Bernard,
CNRS, Laboratoire de Physique, F-69342, Lyon, France

⁶Meteorological Institute, Ludwig Maximilian University Munich, Theresienstr. 37, 80333 Munich, German

Representing the glaciation of mixed-phase clouds in terms of the Wegener-Bergeron-Findeisen process is a challenge for many weather and climate models, which tend to overestimate this process because cloud dynamics and microphysics are not accurately represented. As turbulence is essential for the transport of water vapour from evaporating liquid droplets to ice crystals, we developed a statistical model using established closures to assess the role of small-scale turbulence. The model successfully captures results of direct numerical simulations, and we use it to assess the role of small-scale turbulence. We find that small-scale turbulence broadens the droplet-size distribution somewhat, but it does not significantly affect the glaciation time on submetre scales. However, our analysis indicates that turbulence on larger spatial scales is likely to affect ice growth. While the model must be amended to describe larger scales, the present work facilitates a path forward to understanding the role of turbulence in the Wegener-Bergeron-Findeisen process.

<footnote>* Bernhard.Mehlig@gu.se</footnote>
# I. INTRODUCTION

In mixed-phase clouds, ice particles grow at the cost of evaporating water droplets via the so-called Wegener-Bergeron-Findeisen (WBF) process [1]. This occurs because the saturation vapour pressure differs for liquid water and ice. In the absence of any other processes, the WBF process turns mixed-phase clouds into ice clouds, which is commonly referred to as glaciation. However, mixed-phase clouds can be astonishingly stable (see e.g. [2]), evading a too simplistic interpretation of the WBF process. Thus, our understanding and ability to adequately represent the WBF process has important implications on the longevity and coverage of mixed-phase clouds and hence Earth's radiation budget (e.g. Ref. [3]).

Korolev [4] used a rising-parcel model to study growth of ice particles in mixed-phase clouds, under the assumption that ice particles, water droplets, water vapour and temperature are well mixed, showing that the cooling rates associated with a sufficiently strong updraft can prevent full glaciation. In a similar framework, Ervens et al. [5] highlighted the importance of the number of ice particle on the glaciation process, where smaller ice particle concentrations slow down glaciation.

A question that has received little attention in connection with the WBF process is the impact of small-scale turbulence. Recently, Chen et al. [6] used direct numerical simulations (DNSs) to study the influence of turbulence on growth of ice particles in a mixed-phase cloud, determining which conditions favour ice growth in cloud-top generating cells (CTGC). The authors found that a higher liquid-water content (LWC) and higher relative humidity (RH) favour ice growth by the WBF process. Once the water droplets have evaporated, ice particles continue to grow consuming the remaining water vapour in the cloud. The simulations show that small-scale turbulence has a weak effect on the change of the mean radii of water droplets and ice particles, on the LWC and the ice-particle mass, and therefore on the glaciation time (defined as the time it takes to reach an ice-mass fraction of 0.9). On the other hand, the simulations show how small-scale turbulence increases the width of the particle-size distributions.

Wang et al. [7] compiled parameter sets and specifications to compare models for mixed-phase processes in the Michigan Pi cloud chamber [8], including small-scale turbulence. The corresponding microscopic equations for the core region of the Pi chamber are similar to those of Chen et al. [6], except that a constant influx and removal of ice particles and droplets due to settling is specified.

Here, we analyse a statistical model for these processes, derived from the mapping-closure approximation [9–11] under the assumption that the Lagrangian supersaturation distributions (the distribution of supersaturation along droplet paths) are Gaussian, generalising statistical models for droplet-phase change in turbulence [12–16] to mixed-phase clouds. A strength of the model is that it is constructed from the microscopic governing equations. To assess the accuracy of the model, we compare its predictions to DNS results, for the parameters of Ref. [6], and for the Pi-chamber test case [7]. Another advantage of the model is that it is straightforward to disregard small-scale turbulence in the model, simply by ignoring the stochastic terms. In this way the model simplifies to a parcel model for mixed-phase clouds [4, 17, 18]. This allows us to study under which circumstances small-scale turbulence matters for glaciation, and when it does not. Using the model we investigate the non-dimensional parameters of the problem, and discuss how our conclusions depend on the spatial scale of the turbulent fluctuations.

The remainder of this article is organised as follows. In Section II we describe the microscopic model for mixed-phase clouds, the basis for our DNS, and those of Refs. [6] and [7]. The statistical model is introduced in Section III. Section IV summarises our results, for DNS, statistical model, and for its deterministic limit that disregards small-scale turbulence. In Section V we compare the results and discuss their implications for glaciation of mixed-phase clouds. We summarise our conclusions in Section VI. Three Appendices contain a summary of all parameters used in the calculations, and mathematical details regarding our statistical-model-analysis.

# II. MICROSCOPIC MODEL

## A. Supersaturation over ice and water

The microscopic model of Chen et al. [6] describes how local temperature and the water-vapour mixing ratio are advected by the turbulent flow, and how droplets and ice particles grow and shrink in response to local fluctuations of these quantities. We start by showing how to simplify this microscopic dynamics by combining the water-vapour and temperature fields into supersaturation fields over ice and water. Thereby we extend the results of Refs. [14, 19, 20] that treat the case without ice. Water-vapour supersaturations $s_w$ and $s_i$ over liquid water and ice are defined via the partial pressure $p_v$ of water vapour, and saturated vapour pressures $p_{v,w}$ and $p_{v,i}$ over liquid water and ice:

$$s_w = \frac{p_v}{p_{v,w}} - 1, \quad s_i = \frac{p_v}{p_{v,i}} - 1. \qquad (1)$$
For brevity, we refer to “liquid water” as “water” and use subscripts “w”, “i” and “v” respectively for liquid water, ice and vapor. Expressions for $p_w(T)$ and $p_i(T)$ as functions of temperature $T$ are given by Eq. (A1) in Appendix A. Supersaturations are often expressed in terms of the mixing ratio. The mixing ratio $q_v$ of water vapour is defined as the ratio of the mass $m_v$ of water vapour to the mass $m_a$ of the dry air in a given volume, $q_v = m_v/m_a$; or in terms of densities $q_v = \rho_v/\rho_a$, where $\rho_a = p_a/(R_a T)$ is the dry air density at partial air pressure $p_a$. Note that the full pressure $p$ of the mixture is the sum of partial pressures of air and water vapour, $p = p_v + p_a$. But because $m_v \ll m_a$, we can take $\rho_a \approx \rho$ (the density of the mixture) and $p_a \approx p$, so that $p_v = (R_v/R_a)pq_v$, $R_a$ and $R_v$ being the specific gas constants for dry air and water. Using this relation in Eqs. (1), one can compute the supersaturation over water and ice as

$$s_w = \frac{R_v}{R_a} \frac{p}{p_{v,w}} q_v - 1, \qquad s_i = \frac{R_v}{R_a} \frac{p}{p_{v,i}} q_v - 1. \qquad (2)$$

In order to derive a consistent diffusion-convection-reaction equation for supersaturation, one approximates supersaturation as a linear function of $q_v$, $T$ and $p$ near their reference values $q_{v,0}$, $T_0$ and $p_0$ [14]. To this end we compute the differential of $s_w$ from Eq. (2):

$$ds_w = (1 + s_w) \left( \frac{dq_v}{q_v} - \frac{L_w}{R_v T} \frac{dT}{T} + \frac{dp}{p} \right) = \frac{R_v}{R_a} \frac{p}{p_{v,w}} dq_v + (1 + s_w) \left( -\frac{L_w}{R_v T} \frac{dT}{T} + \frac{dp}{p} \right), \qquad (3)$$

where $L_w$ is the latent heat of water evaporation, and $L_i$ is the latent heat of ice sublimation. We determine $L_w$ and $L_i$ in the following way to ensure consistency with the approximations of $p_{v,w}$ and $p_{v,i}$,

$$L_w(T) = R_v T^2 \frac{d \ln p_{v,w}}{dT}, \qquad L_i(T) = R_v T^2 \frac{d \ln p_{v,i}}{dT}. \qquad (4)$$

We can further simplify Eq. (3). Within the Oberbeck-Boussinesq approximation, the variations of $p$ and $T$ around $p_0$ and $T_0$ are small, allowing us to use constant coefficients in front of the differentials. Next, we deal with the factor $1 + s_w$. We assume that the supersaturation variations $\Delta s_w$ satisfy $\Delta s_w \ll 1 + s_{w,0}$, where

$$s_{w,0} = s_w(q_{v,0}, T_0, p_0) = \frac{R_v}{R_a} \frac{p_0}{p_{v,w}(T_0)} q_{v,0} - 1. \qquad (5)$$

We stress that this assumption is violated when the variations of $s_w$ (equivalently of $q_v$) are not small, as for example when completely dry air at $s_w = -1$ saturates to $s_w = 0$.

When the supersaturation fluctuations are small enough, we can integrate Eq. (3) using the simplifying assumptions of the previous paragraph to obtain $s_w$ as a linear function of $q_v$, $T$ and $p$ (and $s_i$ is derived in a similar manner):

$$s_w = s_{w,0} + \frac{R_v}{R_a} \frac{p_0}{p_{v,w}(T_0)} (q_v - q_{v,0}) + (1 + s_{w,0}) \left( -\frac{L_w(T_0)}{R_v T_0} \frac{T-T_0}{T_0} + \frac{p-p_0}{p_0} \right), \qquad (6a)$$

$$s_i = s_{i,0} + \frac{R_v}{R_a} \frac{p_0}{p_{v,i}(T_0)} (q_v - q_{v,0}) + (1 + s_{i,0}) \left( -\frac{L_i(T_0)}{R_v T_0} \frac{T-T_0}{T_0} + \frac{p-p_0}{p_0} \right). \qquad (6b)$$

In order to derive the diffusion-convection-reaction equations for the fields $s_w(x, t)$ and $s_i(x, t)$ from Eqs. (6), we follow [19] and start from the corresponding equations for the fields $T(x, t)$ and $q_v(x, t)$:

$$\frac{dT}{dt} = \kappa_T \frac{\partial^2 T}{\partial x_j \partial x_j} + \frac{L_w(T_0)}{c_p} C_w + \frac{L_i(T_0)}{c_p} C_i, \qquad (7a)$$

$$\frac{dq_v}{dt} = \kappa_{q_v} \frac{\partial^2 q_v}{\partial x_j \partial x_j} - C_w - C_i. \qquad (7b)$$

Here $C_w(x, t)$ and $C_i(x, t)$ are water and ice condensation and deposition rates, which are discussed in more detail in Section II C, $c_p = \frac{7}{2}R_a$ is the specific heat of air at constant pressure, $\kappa_T$ and $\kappa_{q_v}$ are the molecular diffusivities of $T$ and $q_v$. Math-style Latin indices (as in $x_j$) denote vector/tensor components in Cartesian coordinates, and we use the Einstein summation convention for repeated indices. Next, $d/dt = \partial/\partial t + u_j \partial/\partial x_j$ are the components of the convective derivative, where $u_j(x, t)$ is the turbulent velocity of air determined by the Navier-Stokes equations

$$\frac{du_j}{dt} = -\frac{1}{\rho_0} \frac{\partial p}{\partial x_j} + \nu \frac{\partial^2 u_j}{\partial x_k \partial x_k} + f_j^{(u)}, \qquad \frac{\partial u_j}{\partial x_j} = 0, \qquad (8)$$
Here $\nu$ is the kinematic viscosity of air. The forcing $f_j^{(u)}$ is required to maintain stationary turbulence, as in Refs. [6, 7].

Now we combine Eqs. (6) and (7). For this, we first calculate the Laplacian of $s_w$ from Eq. (6a). Within the Oberbeck-Boussinesq approximation only the hydrostatic pressure affects the thermodynamic variables like $s_w$, and the Laplacian of hydrostatic pressure is negligible compared to the effects of Laplacians of $T$ and $q_v$, which are dominated by small-scale turbulence. Thus, we can neglect the Laplacian $\partial^2 p/\partial x_j \partial x_j$ to obtain

$$ \frac{\partial^2 s_w}{\partial x_j \partial x_j} = \frac{R_v}{R_a} \frac{p_0}{p_{v,w}(T_0)} \frac{\partial^2 q_v}{\partial x_j \partial x_j} - \frac{(1+s_{w,0})L_w(T_0)}{R_v T_0^2} \frac{\partial^2 T}{\partial x_j \partial x_j}. \quad (9) $$

The second step is to assume $\kappa_{q_v} \approx \kappa_T$, which is justified since $\kappa_{q_v} = 1.17\kappa_T$ (Eq. (A11)). Defining $\kappa = \sqrt{\kappa_{q_v} \kappa_T}$ allows us to write

$$ \frac{ds_w}{dt} = \kappa \frac{\partial^2 s_w}{\partial x_j \partial x_j} - A_{2,w,w} C_w - A_{2,w,i} C_i, \quad (10a) $$

$$ \frac{ds_i}{dt} = \kappa \frac{\partial^2 s_i}{\partial x_j \partial x_j} - A_{2,i,w} C_w - A_{2,i,i} C_i, \quad (10b) $$

where the equation for $s_i$ is derived analogously. Here we introduced the parameters

$$ A_{2,\phi_1,\phi_2} = \frac{R_v}{R_a} \frac{p_0}{p_{v,\phi_1}(T_0)} + \frac{(1+s_{\phi_1,0})L_{\phi_1}(T_0)L_{\phi_2}(T_0)}{c_p R_v T_0^2}. \quad (11) $$

In this expression, the variable $\phi$ stands for a particular condensed phase, either $\phi = w$ or $\phi = i$. If we disregard the ice phase, we obtain the supersaturation dynamics used in the statistical models for evaporation of water droplets at the cloud edge [11, 14]. There are minor differences in the expressions of the A-parameters in those references, reflecting slightly different assumptions.

The model (10) can be further simplified if we express $s_i$ as a function of $s_w$, so that it is sufficient to solve a single partial differential equation for $s_w$. The final model for $s_w$ and $s_i$ used in our DNS and in the statistical model is

$$ \frac{ds_w}{dt} = \kappa \frac{\partial^2 s_w}{\partial x_j \partial x_j} - A_{2,w}C_w - A_{2,i}C_i + f^{(s_w)}, \quad (12a) $$

$$ s_i = A_4(s_w + 1) - 1. \quad (12b) $$

Here we introduced a forcing term $f^{(s_w)}$ representing the forcing of small-scale fluctuations due to large spatial scales that are not resolved by the DNS [6, 7]. The new A2-parameters are given by $A_{2,w} = A_{2,w,w}$ and $A_{2,i} = A_{2,w,i}$, in simplified notation. The parameter $A_4$ is defined as:

$$ A_4 = \frac{p_{v,w}(T_0)}{p_{v,i}(T_0)} \quad (13) $$

Since $p_{v,i} > p_{v,w}$ for $T < 0$ °C, we have $A_4 > 1$ and $s_i > s_w$ in the mixed-phase cloud. The WBF process corresponds to $s_i > 0 > s_w$, under this condition water droplets evaporate and ice particles grow. To assess how the single-supersaturation approximation works, we estimate the error $\delta s_i$ between $s_i$ calculated from (6b) and (12b). For $T_0$ between 231.15 K and 273.15 K and the parameters from Appendix A we find:

$$ \delta s_i / s_i \approx -0.011 \text{ K}^{-1} (T - T_0) / s_i. \quad (14) $$

We conclude that the single-supersaturation approximation (12b) works well for temperature variations of the order of $T - T_0 \sim 1$ K, provided that $|s_i| \gg 0.01$. The latter condition is satisfied at the initial and most interesting stage of glaciation, when air is saturated with respect to water, $s_w = 0$ and $s_i > 0.1$. More generally, the approximation works well when temperature fluctuations are much smaller than 1 K, as in the case of the CTGC [6].

## B. Particle dynamics

Water droplets and ice particles are assumed to be so small that they follow the flow, their positions $\mathbf{x}_w(t)$ and $\mathbf{x}_i(t)$ obey

$$ \frac{dx_{w,j}}{dt} = u_j(\mathbf{x}_w, t), \quad \frac{dx_{i,j}}{dt} = u_j(\mathbf{x}_i, t). \quad (15) $$
In other words, effects of particle inertia [21] are neglected. The particle radii $r_w(t)$ and $r_i(t)$ change according to:

$$ \frac{dr_w^2}{dt} = 2A_{3,w} a_3(r_w/r_{A_{3,w}}) [s_w - s_{w,K}(r_w)], \quad (16a) $$

$$ \frac{dr_i^2}{dt} = \begin{cases} 2A_{3,i} a_3(r_i/r_{A_{3,i}}) s_i & \text{if } r_i > 0, \\ 0 & \text{if } r_i = 0, \end{cases} \quad (16b) $$

with supersaturation taken at the particle position, e.g. $s_w(t) = s_w(x_w, t)$. For ice particles, $dr_i/dt$ is constrained to vanish at $r_i = 0$ to ensure that $r_i^2$ remains non-negative, as it must.

Water droplets are not allowed to completely evaporate, due to the Köhler correction term in Eq. (16), involving the radius-dependent function $s_{w,K}$. This function is parameterized by the dry aerosol radius $r_{dry}$ and the hygroscopicity coefficient $\kappa$:

$$ s_{w,K}(r_w) = \frac{r_w^3 - r_{\text{dry}}^3}{r_w^3 - r_{\text{dry}}^3(1 - \kappa)} - 1. \quad (17) $$

A more general expression for Köhler corrections contains an exponential term for Kelvin curvature effects [22]. In Eq. (17) we approximated the exponential by unity. This is a good approximation for our values of $r_{dry}$. Unlike water droplets, ice particles are not allowed to reactivate: once an ice particle evaporates and $r_i$ reaches zero, it stays evaporated with $r_i = 0$ (the ice is pure). This is in agreement with the specifications of both the CTGC and the Pi Chamber cases. Equations (16) contain corrections for the efficiency of accommodation of water vapour on the particle surface, introducing the particle-size dependent function $a_3(x) = x/(x+1)$ with accommodation length $r_{A_{3,\phi}}$ [23, 24]:

$$ r_{A_{3,\phi}} = A_{3,\phi} \frac{\rho_{\phi} \sqrt{2\pi R_v T_0}}{\alpha_{q_v,\phi} p_{v,\phi}(T_0)}, \quad (18) $$

where $\phi$ refers to the phase of condensed water, either $\phi = w$ or $\phi = i$. Particles with radii $r_\phi \gg r_{A_{3,\phi}}$ are not affected by these corrections. Here $\alpha_{q_v,w}$ and $\alpha_{q_v,i}$ are water-vapour accommodation coefficients over water and ice. The use of the particular forms of $A_{3,\phi}$, $a_3$ and $r_{A_{3,\phi}}$ for our cases is justified in Appendix A. Neglecting all radius-dependent corrections corresponds to $a_3 = 1$ and $s_{w,K} = 0$. The A$_3$-parameters in Eqs. (16) and (18) are given by [25]

$$ A_{3,\phi} = \left[ \frac{R_a \rho_{\phi} L_{\phi}^2(T_0)}{R_v \kappa_T c_p T_0 p_0} + \frac{\rho_{\phi} R_v T_0}{\kappa_{q_v} p_{v,\phi}(T_0)} \right]^{-1}. \quad (19) $$

Finally, the specifications for the cloud-chamber test case allow for injection and removal of water droplets and ice particles [7]. Water droplets are injected in the form of dry aerosol at a constant rate $I_w [(\text{m}^3)^{-1}]$ and removed at a rate defined by the settling velocity $u_{\infty,w}$. Ice particles are treated the same way, with $I_i$ and $u_{\infty,i}$. Namely, each particle is removed with a probability $P_w$ (droplet) or $P_i$ (ice particle):

$$ P_w = \min\left(\frac{u_{\infty,w} \Delta t}{H}, 1\right), \quad P_i = \min\left(\frac{u_{\infty,i} \Delta t}{H}, 1\right) \quad (20) $$

where $H$ is the total height of the Pi chamber. The settling velocities depend on the particle radii $r_w$ and $r_i$:

$$ u_{\infty,w} = k_{\infty,w} r_w^2, \quad u_{\infty,i} = k_{\infty,i} r_i^2. \quad (21) $$

The values of the parameters $k_{\infty,w}$ and $k_{\infty,i}$ are specified by [7]. We note that particle shape impacts the sedimentation velocity. Here we assume spherical particles, but larger ice crystals (with radii $> 30\mu\text{m}$) tend to be non-spherical. This is not accounted for in the model (for the data discussed below, Figures 1 and 2, the ice particles do not exceed this size). The overall number of particles changes as

$$ \frac{dN_w}{dt} = V I_w + \frac{N_w}{H} \langle u_{w,\infty} \rangle, \quad \frac{dN_i}{dt} = V I_i + \frac{N_i}{H} \langle u_{i,\infty} \rangle, \quad (22) $$

where $V$ is the volume of the simulation domain (corresponding to the core of the Pi chamber).
### C. Condensation and deposition rates

The condensation and deposition rates $C_w$, $C_i$ for water and ice are defined through the rate of change of condensed water content:

$$C_w(\mathbf{x}, t) = \frac{4}{3}\pi \frac{\rho_w}{\rho_0} \sum_{\alpha=1}^{N_w} G(\mathbf{x} - \mathbf{x}_{w,\alpha}) \frac{dr_{w,\alpha}^3}{dt}, \qquad (23a)$$

$$C_i(\mathbf{x}, t) = \frac{4}{3}\pi \frac{\rho_i}{\rho_0} \sum_{\alpha=1}^{N_i} G(\mathbf{x} - \mathbf{x}_{i,\alpha}) \frac{dr_{i,\alpha}^3}{dt}, \qquad (23b)$$

where $N_w$ is the number of water droplets, $N_i$ is the number of ice particles, and $G$ is the standard spatial kernel, normalized to unity [6, 26]. The spatial range of $G$ is the linear size of a DNS-grid cell. Using Eqs. (16) we can rewrite $C_w$ and $C_i$ as

$$C_w(\mathbf{x}, t) = 4\pi \frac{\rho_w}{\rho_0} A_{3,w} \sum_{\alpha=1}^{N_w} G(\mathbf{x} - \mathbf{x}_{w,\alpha}) r_{w,\alpha} a_3(\frac{r_{w,\alpha}}{r_{A_{3,w}}}) [s_{w,\alpha} - s_{w,K}(r_{w,\alpha})], \qquad (24a)$$

$$C_i(\mathbf{x}, t) = 4\pi \frac{\rho_i}{\rho_0} A_{3,i} \sum_{\alpha=1}^{N_i} G(\mathbf{x} - \mathbf{x}_{i,\alpha}) r_{i,\alpha} a_3(\frac{r_{i,\alpha}}{r_{A_{3,i}}}) s_{i,\alpha}. \qquad (24b)$$

Here $s_{w,\alpha}$ is the supersaturation field at the position of particle $\alpha$, $s_{w,\alpha}(t) = s_w(\mathbf{x}_\alpha, t)$. Note that the multiplication by $r_{i,\alpha}$ in Eq. (24b) correctly accounts for the $r_i = 0$ condition in Eq. (16b).

Without injection and sedimentation ($N_w$ and $N_i$ are constant) and for vanishing mean forcing $\langle f^{(s_w)} \rangle$, Eqs. (12a) and (23) imply that the quantities

$$s_{w,inv} = \langle s_w \rangle_V + \frac{4}{3}\pi \frac{\rho_w}{\rho_0} A_{2,w} n_w \langle r_w^3 \rangle_V + \frac{4}{3}\pi \frac{\rho_i}{\rho_0} A_{2,i} n_i \langle r_i^3 \rangle_V, \qquad (25a)$$

$$s_{i,inv} = A_4(s_{w,inv} + 1) - 1. \qquad (25b)$$

are invariant. Here $n_w = N_w/V$ and $n_i = N_i/V$ are droplet and ice number densities, and $\langle s_w \rangle_V = \frac{1}{V} \int s_w d\mathbf{x}$ is the spatially averaged supersaturation. Since Eq. (12a) for $s_w$ is derived from the two Eqs. (7) for $q_v$ and for $T$, the conservation law (25) combines both water and thermal energy conservation during condensation and deposition of water vapor [13, 14]. Physically, $s_{w,inv}$ and $s_{i,inv}$ correspond to supersaturation over water and ice if all the particles evaporate, so that all the water is contained in the form of water vapor. In case we ignore the Köhler corrections, $s_{i,inv}$ provides us with an insight into the final state of the cloud. Since water evaporates ($r_w = 0$), $s_{i,inv} > 0$ means that ice remains, while $s_{i,inv} \le 0$ implies that ice evaporates too.

All thermodynamic parameters in the above microscopic equations are summarised in Tables IV (CTGC) and V (Pi chamber) in Appendix A.

### D. Direct numerical simulations

We performed DNSs using Eqs. (8) to (22) for mixed-phase processes in the core of the Pi chamber, for the parameters specified by Wang et al. [7]. The turbulent dissipation rate per unit mass was $\epsilon = 66 \text{ cm}^2 \text{s}^{-3}$ in a cubic domain with side length $L = 20$ cm. With $\nu = 1.278 \times 10^{-5} \text{ m}^2 \text{s}^{-1}$, the Kolmogorov length $\eta = (\nu^3/\epsilon)^{1/4}$ is around $\eta = 0.75$ mm. To properly resolve the turbulent flow we used a numerical resolution with 256³ collocation points to solve the Eulerian equations.

As specified by Wang et al. [7], our DNSs used the forcing term in the Navier-Stokes equations (8) to maintain a statistically steady turbulent state with constant dissipation rate $\epsilon$. In Fourier space, the forcing reads

$$\hat{f}^{(u)}(\mathbf{k}) = \epsilon N(t) \hat{u}(\mathbf{k}, t) \quad \text{for } |\mathbf{k}| < 3\pi/L. \qquad (26)$$

Here $\hat{u}(\mathbf{k}, t)$ is the Fourier transform of the turbulent velocity field $u(\mathbf{x}, t)$, $\mathbf{k}$ is the wave vector, and $N(t) = (\sum_{|\mathbf{k}|<2\pi/L} |\hat{u}(\mathbf{k}, t)|^2)^{-1}$ is a normalisation factor. The supersaturation equation (12a) is also forced. For the CTGC we use a Gaussian random forcing [27],

$$\hat{f}^{(s_w)}(\mathbf{k}) = \beta \text{dW}(\mathbf{k}, t) \quad \text{for } |\mathbf{k}| < 3\pi/L, \qquad (27)$$
<table><caption>TABLE I. DNS time- and length scales for the Pi chamber [7].</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Simulation time t<sub>DNS</sub></td><td>600 s</td></tr><tr><td>Eddy-turnover time k/ε</td><td>1.6 s</td></tr><tr><td>Kolmogorov time τ<sub>η</sub></td><td>4.4 × 10<sup>-2</sup> s</td></tr><tr><td>Integration time step Δt</td><td>1 × 10<sup>-3</sup> s</td></tr><tr><td>Linear domain size L</td><td>2 × 10<sup>-1</sup> m</td></tr><tr><td>Integral length scale L<sub>int</sub></td><td>4.2 × 10<sup>-2</sup> m</td></tr><tr><td>Taylor microscale λ</td><td>1.6 × 10<sup>-2</sup> m</td></tr><tr><td>Kolmogorov length η</td><td>7.5 × 10<sup>-4</sup> m</td></tr><tr><td>Spatial resolution Δx</td><td>7.8 × 10<sup>-4</sup> m</td></tr></tbody></table>

where $dW(t)$ is white noise with unit variance independently chosen for different $k$. The numerical factor $\beta$ is used to maintain the prescribed steady water supersaturation root mean square $\sigma_{s_w}$ before the aerosol injection. The forcing (27) differs from the one suggested by Wang et al. [7] where there is no explicit forcing term, instead the supersaturation Fourier coefficients of the forced wavenumbers are rescaled at each time step to maintain the prescribed $\sigma_{s_w}$. We tested that the two forcing schemes yield the same condensation/evaporation statistics, provided that they achieve the same statistically steady-state value of $\sigma_{s_w}$. For the Pi chamber, Wang et al. [7] specify an additional average forcing in Eq. (12a) that nudges the average supersaturation,

$$ \langle f^{(s_w)} \rangle = - \left( \langle s_w \rangle - s_{w,\text{force}} \right) / \tau_{s_w,\text{force}}, \quad (28) $$

where $s_{w,force}$ is the mean supersaturation before aerosol injection, and $\tau_{s_w,force}$ is a forcing timescale. This forcing mimics the property of the Rayleigh–Benard convection inside the cloud chamber to achieve a statistically steady thermodynamic state [28].

For the Pi chamber, we need to add and remove the particles during each time step $\Delta t$ as specified by Wang et al. [7], see Eqs. (21) and (22) in Section II B. To this end, the numbers of added particles with radii $r_{w,initial}$ and $r_{i,initial}$ are

$$ \Delta N_w = \text{floor}(I_w V \Delta t) + \begin{cases} 1, & \text{frac}(I_w V \Delta t) \ge \xi_w, \\ 0, & \text{else,} \end{cases} \quad (29) $$

$$ \Delta N_i = \text{floor}(I_i V \Delta t) + \begin{cases} 1, & \text{frac}(I_i V \Delta t) \ge \xi_i, \\ 0, & \text{else,} \end{cases} \quad (30) $$

where $I_w$ and $I_i$ are water and ice injection rates, and $\xi_w$ and $\xi_i$ are independent random variables uniformly distributed in [0,1]. The initial ice particle radius is 2 $\mu$m, while the initial droplet size correspond to a dry aerosol particle with a diameter of 0.125 $\mu$m.

For our DNSs, we used the same numerical solver as in Refs. [14, 19]. The Navier-Stokes equations (8) were solved in Fourier space using fast Fourier transform. The nonlinear terms were calculated in configuration space using the de-aliasing 2/3 rule [19]. Time integration used a low-storage third-order Runge-Kutta method, where the terms are treated exactly by using integration factors, while the non-linear terms followed an Adam-Bashforth scheme. The same Runge-Kutta scheme was used to integrate the equations of motion (15) for water droplets and ice particles, and their growth equations (16). A linear interpolation scheme was used to evaluate the air velocity and supersaturation at the particle positions, while linear extrapolation was employed to calculate the condensation rates $C_w$ and $C_i$ in Eqs. (24). The parameter values for the DNSs are summarised in Table I.

We also performed our own DNSs for the CTGC [6], for the same parameters as in Ref. [6]. They advect two scalar fields, temperature and water-vapour mixing ratio (Section V). The parameters for these DNS runs are given in Table II. We note that we used slightly larger times steps and slightly coarser grid than Chen et al. [6]. We use our DNSs for the CTGC to determine the Lagrangian correlation time of supersaturation, an input needed for the statistical model that is discussed next.

### III. STATISTICAL MODEL

To understand glaciation dynamics and how it is affected by small-scale turbulence one could simulate the mi- croscopic model described in Section II for a wide range of parameters. Here we take an alternative approach: we derive a statistical model that allows us to systematically study the parameter dependencies of the glaciation process,
<table><caption>TABLE II. DNS time- and length scales for the CTGC [6].</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Simulation time t<sub>DNS</sub></td><td>95 s</td></tr><tr><td>Eddy turnover time k/ε</td><td>2.82 s</td></tr><tr><td>Kolmogorov time τ<sub>η</sub></td><td>1.26 × 10<sup>-1</sup> s</td></tr><tr><td>Integration time step Δt</td><td>2.5 × 10<sup>-3</sup> s</td></tr><tr><td>Linear domain size L</td><td>2 × 10<sup>-1</sup> m</td></tr><tr><td>Integral length scale L<sub>int</sub></td><td>1.31 × 10<sup>-1</sup> m</td></tr><tr><td>Taylor microscale λ</td><td>2.13 × 10<sup>-2</sup> m</td></tr><tr><td>Kolmogorov length η</td><td>1.42 × 10<sup>-3</sup> m</td></tr><tr><td>Spatial resolution Δx</td><td>1.5625 × 10<sup>-3</sup> m</td></tr></tbody></table>

and provides immediate insight into possible effects of small-scale turbulence. We validate the model by showing that it yields quantitative agreement with the DNS results of [6]. In essence, the model is a statistical model for the supersaturation, approximating Eqs. (8, 12, 15), while the Eqs. (16) for particle radii remain the same. The derivation of the model rests on two assumptions:

A1. The supersaturation statistics along water-droplet, ice-particle, and Lagrangian fluid paths are the same.

A2. The supersaturation statistics are Gaussian.

Our DNSs for the CTGC and for the core of the Pi chamber show that these assumptions hold. In Section V we discuss their range of validity. To derive the model under the above assumptions, we start by decomposing the supersaturation along a particle trajectory into its mean and fluctuating parts

$$s_w = \langle s_w \rangle + s'_w. \qquad (31)$$

We use the usual notation $\langle \cdot \rangle$ for ensemble averages of physical quantities and $\cdot'$ for their fluctuating parts. The system is statistically homogeneous: the mean values may depend on $t$, but not on $x$. Since water droplets and ice particles are Lagrangian tracers, and since the flow is incompressible, single-point Eulerian and Lagrangian statistics are the same. Therefore we can take the ensemble average of Eq. (12a) to obtain the evolution equation for the mean supersaturation $\langle s_w \rangle$:

$$\frac{d\langle s_w \rangle}{dt} = -A_{2,w}\langle C_w \rangle - A_{2,i}\langle C_i \rangle + \langle f^{(s_w)} \rangle. \qquad (32)$$

To close Eq. (32), we need expressions for the mean condensation rates $\langle C_w \rangle$ and $\langle C_i \rangle$, which we derive in Appendix B:

$$\langle C_w \rangle = \frac{4}{3} \pi \frac{\rho_w}{\rho_0} n_w \left\langle \frac{dr_w^3}{dt} \right\rangle, \quad \langle C_i \rangle = \frac{4}{3} \pi \frac{\rho_i}{\rho_0} n_i \left\langle \frac{dr_i^3}{dt} \right\rangle. \qquad (33)$$

Using Eqs. (16), these expressions evaluate to

$$\langle C_w \rangle = 4\pi \frac{\rho_w}{\rho_0} A_{3,w} n_w \langle a_3 \left( \frac{r_w}{r_{A_3,w}} \right) r_w \left[ (s_w - s_{w,K}(r_w)) \right] \rangle, \qquad (34a)$$

$$\langle C_i \rangle = 4\pi \frac{\rho_i}{\rho_0} A_{3,i} n_i \langle a_3 \left( \frac{r_i}{r_{A_3,i}} \right) r_i s_i \rangle. \qquad (34b)$$

Note that averages involving particle radii also include averaging over particles, e.g. $\langle r_w s_w \rangle = \frac{1}{N_w} \sum_{\alpha=1}^{N_w} \langle r_{w,\alpha} s_{w,\alpha} \rangle$. However for simplicity we do not introduce a special notation for such averages, except for Appendix B where we use $\langle \cdot \rangle_w$ or $\langle \cdot \rangle_i$. For the Pi chamber, the average forcing $\langle f^{(s_w)} \rangle$ in (32) is given by Eq. (28). For the CTGC [6], the average vanishes.

With a model for the mean supersaturation in place, we now introduce a model for its fluctuating part. Fries et al. [11] used the mapping closure of Pope [9] and Chen et al. [10] to accurately describe non-Gaussian dynamics of $s'_w$ during the evaporation of water droplets at the cloud edge. We start from the same model here. Since in our case the statistics of $s_w$ is Gaussian, the mapping closure becomes a linear theory and reduces to an Ornstein–Uhlenbeck process for $s'_w$

$$ds'_{w} = -\frac{1}{\tau_{s_w}^{(L)}} s'_{w} dt + \sqrt{\frac{2\sigma_{s_w}^2}{\tau_{s_w}^{(L)}}} dW(t). \qquad (35)$$
Here $dW(t)$ are white-noise increments, while the supersaturation variance $\sigma_{s_w}^2 = \langle s_w'^2 \rangle$ and the correlation time $\tau_{s_w}^{(L)}$ are the two parameters of the model. The model (35) is also known as the Langevin mixing model [9, Eq. 5.52]. Related models have been used to describe the effect of supersaturation fluctuations on the growth of water droplets in turbulent clouds [12, 13, 15, 16, 29]. Some of them contain additional condensation terms in the equation for $s_w'$. To understand why such terms do not matter in our case, consider a more general statistical model for the supersaturation fluctuations

$$ds_w' = -A_{2,w} \langle C_w' | s_w', t \rangle - A_{2,i} \langle C_i' | s_w', t \rangle - \frac{s_w'}{\tau_{s_w}^{(L)}} dt + \sqrt{D^{(2)}} dW(t). \quad (36)$$

Here $\langle C'_\phi | s'_w, t \rangle = \langle C_\phi | s'_w, t \rangle - \langle C_\phi \rangle$, where $\langle C_\phi | s_w, t \rangle$ are conditional condensation rates, and $D^{(2)}$ is chosen such that $\sigma_{s_w}^{(2)}$ [Eq. (B16)]. Appendix B outlines how to derive Eq. (36)) using the method of Sarnitsky and Heinz [30].

However, the fluctuating condensation-rate contributions are negligible if the timescale of turbulent mixing $\tau_{s_w}^{(L)}$ is much smaller than the timescales $\tau_{s_w,w}$ and $\tau_{s_w,i}$ of supersaturation evolution due to the phase change of water droplets and ice particles,

$$\tau_{s_w,w} = \frac{\sigma_{s_w}^2}{A_{2,w} |\langle C_w' s_w' \rangle|}, \quad \tau_{s_w,i} = \frac{\sigma_{s_w}^2}{A_{2,i} |\langle C_i' s_w' \rangle|}, \quad (37)$$

where $|\cdot|$ denotes the absolute value. Their fractions with $\tau_{s_w}^{(L)}$ define the supersaturation Damköhler numbers

$$Da_{s_w,w} = \frac{\tau_{s_w}^{(L)}}{\tau_{s_w,w}}, \quad Da_{s_w,w} = \frac{\tau_{s_w}^{(L)}}{\tau_{s_w,w}}. \quad (38)$$

We conclude: for small Damköhler numbers, $Da_{s_w,w} \ll 1$ and $Da_{s_w,i} \ll 1$, one can use the model (35) instead of (36). Our model calculations confirm that the supersaturation Damköhler numbers are smaller than unity for the cases studied here.

To close the model (35), we need to provide the supersaturation variance $\sigma_{s_w}^2$ and the correlation time $\tau_{s_w}^{(L)}$. For the Pi chamber, the variance is given by Wang et al. [7]. For the CTGC, Chen et al. [6] specify the variances of temperature and the water-vapour mixing ratio, which allows us to compute $\sigma_{s_w}^2$. The Lagrangian correlation time $\tau_{s_w}^{(L)}$ is defined as

$$\tau_{s_w}^{(L)} = \sigma_{s_w}^{-2} \int_0^\infty dt \langle s'(t)s'(0) \rangle, \quad (39)$$

where $\langle s'(t)s'(0) \rangle$ is the Lagrangian autocovariance, i.e. it is taken along Lagrangian trajectories. The numerical values of $\tau_{s_w}^{(L)}$ and $\tau_L$, are given in Appendix A. We note that for the cases studied here, phase change does not affect $\tau_{s_w}^{(L)}$ since the condensation terms are negligible in the dynamics of $s_w'$ for small supersaturation Damköhler numbers. Disregarding phase change, $\tau_{s_w}^{(L)}$ is commonly related to the large eddy turbulent timescale $\tau_L = k/\epsilon$:

$$\tau_{s_w}^{(L)} = \frac{C_{0,s_w}}{C_{s_w}} \tau_L, \quad (40)$$

The Lagrangian Obukhov–Corrsin constant $C_{0,s_w}$ comes from the Kolmogorov hypothesis extended to passive scalars and connects the correlation timescale $\tau_{s_w}^{(L)}$ to the dissipation timescale $\sigma_{s_w}^2/\epsilon_{s_w}$, $\tau_{s_w}^{(L)} = C_{0,s_w}\sigma_{s_w}^2/\epsilon_{s_w}$, where $\epsilon_{s_w} = 2\kappa\langle \frac{\partial s_w}{\partial x_j} \frac{\partial s_w}{\partial x_j} \rangle$ is the supersaturation dissipation rate. For the CTGC and Pi chamber, the values of $C_{0,s_w}$ are 1.2 and 1.1, inferred from the DNS described in Section II D. The quantity $C_{s_w}$ is the so called mechanical-to-scalar timescale ratio, formally defined as $C_{s_w} = k\epsilon_{s_w}/(\sigma_{s_w}^2\epsilon)$. The quantity $C_{s_w}$ can be considered approximately constant only in specific types of flows, like the forced isotropic turbulence we deal with here [31]. Its numerical values calculated from the DNS for the CTGC and Pi Chamber cases are 1.8 and 2.2 respectively. We stress again this discussion is valid only for $Da_{s_w} \ll 1$ and $Da_{s_i} \ll 1$. In case of non-negligible supersaturation Damköhler numbers, Fries et al. [11] found that both $C_{0,s_w}$ and $C_{s_w}$ (denoted there as $2/C$ and $2\phi_*$) cannot be considered constant.

To numerically integrate the statistical model, we use the Euler–Maruyama scheme with a time step of 0.05 s for the CTGC, and 0.02 s for the Pi chamber. To ensure that the numerical integration conserves the invariants (25), condensation rates are computed directly from Eqs. (33) and not Eqs. (34). For the CTGC, we use $N_w = N_i = 10^7$ to suppress the statistical noise for cases 1 and 2. The choice of $N_w$ and $N_i$ here has no other consequences, since the model depends only on $n_w$ and $n_i$ in Eq. (33), which are fixed in each run. For the Pi chamber, the number of
<figure><img src="image_2.png" /><figcaption>FIG. 1. Model results for ice growth in cloud-top generating cells [6]. Shown are the DNS results of Chen et al. [6] (solid lines) for the mean droplet radius (a), liquid-water content LWC (b), the mean ice-particle radius (c) and the ice mass (d) as functions of time. In each panel, curves for four parameter sets are shown, these parameter sets are given in Table III. Only ice particles with $r_i > 0.001$ µm are included in the statistics [6]. Also shown are simulations of the statistical model (dashed lines), and of its deterministic limit (dash-dotted lines). The deterministic limit is so close to the full statistical-model results that the lines are hard to distinguish.</figcaption></figure>

<table><caption>TABLE III. Parameters for the CTGC case (Fig. 1).</caption><thead><tr><th>Case</th><th>s<sub>i,inv</sub></th><th>initial &lang;s<sub>w</sub>&rang;</th><th>initial r<sub>w</sub>, µm</th><th>initial r<sub>i</sub>, µm</th><th>n<sub>w</sub>, cm<sup>-3</sup></th><th>n<sub>i</sub>, cm<sup>-3</sup></th><th>σ<sub>s<sub>w</sub></sub></th></tr></thead><tbody><tr><td>1</td><td>-0.080</td><td>-0.2</td><td>10</td><td>1</td><td>1</td><td>100</td><td>0.017</td></tr><tr><td>2</td><td>-0.023</td><td>-0.15</td><td>10</td><td>1</td><td>1</td><td>100</td><td>0.017</td></tr><tr><td>3</td><td>0.42</td><td>-0.1</td><td>10</td><td>1</td><td>100</td><td>10</td><td>0.016</td></tr><tr><td>4</td><td>0.42</td><td>-0.1</td><td>10</td><td>1</td><td>100</td><td>100</td><td>0.016</td></tr></tbody></table>

particles is determined by the particle injection and removal process, which is implemented as described in Section II D.

Below, we refer to the deterministic limit of the statistical model, or deterministic model. It is obtained by taking the limit $\sigma_{s_w} \to 0$ in the statistical model, which amounts to removing the white-noise term in Eq. (35) and setting $s'_w = 0$.

IV. RESULTS

Figure 1 shows the results for ice growth and water evaporation due to the WBF in the CTGC [6]. We chose the four cases from the accompanying Replication Data [32] where the particle-size evolution is most rapid, their parameters are listed in Tables III and V in Appendix A. Panel (a) shows how the droplet radius shrinks because the droplets evaporate. The radius saturates at a small value determined by the interplay between solute and curvature effects on the one hand, and evaporation on the other hand [33]. Shown are the DNS results of Chen et al. [6] with two scalar fields, temperature and water-vapour mixing ratio (solid lines). We see that the statistical model results (dashed lines) agree very well with those of the DNS, although the droplets evaporate somewhat faster in the statistical model. Also shown are results for the deterministic limit of the model (dash-dotted lines). They are almost indistinguishable
<figure><img src="image_3.png" /><figcaption>FIG. 2. Model results for ice growth in the core of the Pi chamber [7]. Shown are the DNS results (Section IID) (solid lines) for the mean droplet radius (a), the liquid-water content LWC (b), the mean ice-particle radius (c), the ice mass (d), water droplet concentration (e) and ice particle concentration (f) as functions of time. In each panel, curves for five different ice-particle injection rates [cm⁻³ min⁻¹] are shown, the parameter values are given in the insets. Also shown are simulations of the statistical model (dashed lines), and of its deterministic limit (dash-dotted lines). In most but not all cases, the deterministic limit is so close to the full statistical-model results that the lines are hard to distinguish.</figcaption></figure>

from the full statistical-model results. This shows that turbulence has no effect on the evolution of the mean droplet radius. Panel (c) shows how the LWC decreases as the droplets evaporate, with analogous conclusions. In panel (b), we compare how the mean ice-particle radius changes as a function of time. In agreement with the values of the invariant $s_{i,inv}$ (Eq. (25b), Table III): the first two cases have $s_{i,inv} < 0$ and the ice evaporates, while for the last two cases $s_{i,inv} > 0$ and the ice particles grow. Panel (d) shows the IWC, which approaches a non-zero steady state when the cloud glaciates, but tends to zero in the other two cases, as expected. The fact that LWC decreases and IWC increases for cases 3 and 4 indicates that the ice particles grow at the expense of water droplets, as described by the WBF process. In summary, the main conclusion from Fig. 1 is that small-scale turbulence does not affect the mean particle radius. Chen et al. [6] came to the same conclusion for their base case. Comparing the deterministic limit of our statistical model to their DNS for all parameter settings listed in the replication data [32] show that turbulence does not affect the mean particle radii for any of the cases.

Figure 2 shows mean particle radii versus time for the five Pi chamber cases specified by Wang et al. [7] corresponding
to different ice-particle injection rates. As specified by Wang *et al.* [7], the mean droplet radius in Fig. 1 is computed excluding droplets with radii < 3.5 µm. Panel (a) reveals how the droplet radius changes. For the two highest ice-particle injection rates, the radius tends to zero. In other words, the cloud glaciates. For lower ice-particle injection rates, droplets remain in the centre of the Pi chamber at the end of the simulation. The transition to glaciation occurs for ice-injection rates between 5 and 10 cm⁻³ min⁻¹. As for the CTGC, we see that the statistical model (dashed lines) describes the DNS results (solid lines) very well, as does the deterministic limit. Also here we conclude that small-scale turbulence has little effect, at least on droplets of radii > 3.5 µm. Panel (c) shows how the LWC changes as a function of time. The results indicate that turbulence may change the location of the glaciation transition, as evident from the case *I*ᵢ = 5 cm⁻³ min⁻¹ (green lines). The LWC with turbulence (the dashed line) increases at large times and the cloud remains mixed phase, whereas LWC without turbulence (the dash-dotted line) decreases and the cloud glaciates. Approaching the glaciation transition, LWC curves show a dip, as demonstrated the clearest by the case *I*ᵢ = 3 cm⁻³ min⁻¹ (the orange curve) near *t* ~ 180 s. This is explained as follows. Initially the droplets evaporate, so the LWC falls. But as aerosol is added, the competition for water vapour increases, which leads to smaller droplets. Because these droplets do not sediment substantially, the LWC grows. Hence as the overall number of droplets keeps increasing, the mean radius of droplets decreases and reaches a steady state, as demonstrated by panel (e). Hence injecting ice particles results in a decrease of droplets size, but in an increase of droplet concentration. Panels (c), (d) and (f) show how the mean ice-particle radius approaches a plateau, as does the IWC and the ice particle concentration. The differences in water contents and particle concentrations between the statistical model and the DNS are due to slight differences in the particle injection rates. The general message from Fig. 2 is that small-scale turbulence has little effect on the mean radius of ice particles and droplets with r<sub>w</sub> > 3.5 µm for the parameters from Wang *et al.* [7], except possibly upon the timing of the glaciation transition.

Figure 3 illustrates how the fluctuations in droplet radii develop as a function of time for the CTGC. For all four cases, the relative dispersion $\sigma_{r_w}/\langle r_w \rangle$ ($\sigma_{r_w} = \sqrt{\langle r_w^2 \rangle}$) increases rapidly, before decaying to a plateau (for the green curve the decay is not shown).

We see that the statistical model (dashed lines) describes the DNS (solid lines) very well; the deviations are consistent with the attenuation of $\sigma_{s_w}$ with time in the DNS of [6], while in the statistical model we use $\sigma_{s_w}$ constant in time. More importantly, in the deterministic limit $\sigma_{r_w} = 0$ (not shown). We conclude: the particle-size dispersion is a consequence of small-scale turbulence, and it is well described by the statistical model. We do not show corresponding results for ice, because ice evaporates quickly and $\langle r_i \rangle$ reaches zero for cases 1 and 2. For cases 3 and 4, ice grows so rapidly that the relative fluctuations in the ice-particle radii are negligible.

Figure 4 summarises the same for the Pi chamber, regarding the relative dispersion of droplet radii as a function of time. A major difference in this case is that particle removal and injection causes a particle-size dispersion. This is well described by the deterministic model. So here small-scale turbulence is less important for the particle-size dispersion, compared with the CTGC results in Fig. 3. Note that this observation is valid only for larger droplets with $r_w > 3.5$ µm, this cutoff also being the reason why $\sigma_{r_w} = 0$ for glaciated clouds. Now we turn to the question of how turbulence affects smaller droplets (haze).

Figure 5 shows the probability density functions (PDFs) of particle radii for the Pi chamber at large times, for water droplets [panel (a)] and for ice particles [panel (b)]. For the ice particles, statistical-model predictions and the deterministic limit agree very well, for the droplets also, but not near or after the glaciation transition. The smaller the mean particle radius the greater the shift of the deterministic limit PDF to smaller radii compared to the statistical-model PDF. Thus, we see that turbulence widens the distribution of droplet radii for smaller droplets.

V. DISCUSSION

Figures 1 and 2 show that the statistical-model predictions agree very well with the DNS results. We now explain why this is the case here, and under which circumstances the model may fail. To this end we study the limits of validity of the statistical model, which we derived from the assumptions A1 and A2 listed in Section III.

The first assumption, A1, is that the supersaturation statistics along water-droplet, ice-particle, and Lagrangian fluid paths are the same. We expect this to hold when the Damköhler numbers $Da_{s,w}$ and $Da_{s,i}$ [Eq. (38)] are small, and when water droplets and ice particles are initially well-mixed. Under these conditions, turbulence transports water vapour to and from particles faster than phase change occurs. Hence the particles have no time to form a supersaturation field in their vicinity, which might differ from particle-free regions of the flow; all the particles experience the same supersaturation statistics.

The second assumption, A2, is that the supersaturation statistics are Gaussian. In a homogeneous system (no scalar gradients [34–36]), the steady-state distribution of a passive scalar in isotropic homogeneous turbulence is Gaussian [37–40]. Non-Gaussian tails that may be present in transient mixing [41] disappear as the steady state is approached.

In summary, the statistical model from Section III can be justified when the Damköhler numbers $Da_{s,w}$ and $Da_{s,i}$
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Relative dispersion of droplet radii for ice growth in cloud-top generating cells [6]. Shown are the DNS results of Chen et al. [6] (solid lines), simulations of the statistical model (dashed lines). (b) Statistical-model results for the Damköhler number $Da_{rw}$ versus time for the same cases as shown in panel (a).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Relative dispersion of droplet radii for the Pi chamber [7], for different ice-particle injection rates [cm$^{-3}$ min$^{-1}$]. Shown are the DNS results (Section II D, solid lines), simulations of the statistical model (dashed lines), and of its deterministic limit (dash-dotted lines).</figcaption></figure>

are small. The cases summarised in Section IV all have Damköhler numbers smaller than unity, explaining why the statistical model works so well. On the other hand, it is worth noting that the interaction between turbulence and phase change can be more intricate and harder to describe at larger Damköhler numbers. This is in line with the findings of Fries et al. [11], and requires more refined mapping-closure approximations [9, 10].

Next, we quantify the role of turbulence on the dynamics of mean supersaturation and upon the mean particle radii. In the framework of the statistical model, the question is how $\sigma_{sw}$ – the measure of turbulence intensity – enters the dynamical equations for $\langle s_w \rangle$ and $\langle r_w \rangle$. To this end we need to introduce two more Damköhler numbers [besides the
<figure><img src="image_6.png" /><figcaption>FIG. 5. Final probability distributions of particle radii for the Pi chamber [7], for different ice-particle injection rates [cm$^{-3}$ min$^{-1}$]. Shown are statistical-model results (solid lines), and from its deterministic limit (dashed). (a) water droplets, (b) ice particles.</figcaption></figure>

supersaturation Damköhler numbers (38)]:

$$ \mathrm{Da}_{r_w} = \frac{\tau_{s_w}^{(L)}}{\tau_{r_w}}, \qquad \mathrm{Da}_{r_i} = \frac{\tau_{s_w}^{(L)}}{\tau_{r_i}}. \qquad (41) $$

They are associated with the timescales $\tau_{r_w}$ and $\tau_{r_i}$ of the particle-size evolution:

$$ \tau_{r_w} = \frac{\langle r_w \rangle^2}{\left| \frac{d}{dt} \langle r_w^2 \rangle \right|}, \qquad \tau_{r_i} = \frac{\langle r_i \rangle^2}{\left| \frac{d}{dt} \langle r_i^2 \rangle \right|} \qquad (42) $$

We note that these Damköhler numbers scale as $\sim 2A_3s_\phi\tau_{s_w}^{(L)}/\langle r_\phi \rangle^2$ ($\phi = w$ or $\phi = i$). Therefore they are larger for smaller particles.

Figures 1 and 2 imply that small-scale turbulence has only a weak effect on the evolution of the average supersaturation. To explain this, we look at the evolution equation (32) for $\langle s_w \rangle$. It is written in terms of the mean condensation rates $\langle C_w \rangle$ and $\langle C_i \rangle$. For the CTGC (no particle injection or removal), we can estimate $\langle C_w \rangle$ and $\langle C_i \rangle$ as

$$ \langle C_{\phi} \rangle = 4\pi \frac{\rho_{\phi}}{\rho_0} A_{3,\phi} n_{\phi} (\langle r_{\phi} \rangle \langle s_{\phi} \rangle + \langle r'_{\phi} s'_{\phi} \rangle) \approx 4\pi \frac{\rho_{\phi}}{\rho_0} A_{3,\phi} n_{\phi} \langle r_{\phi} \rangle \langle s_{\phi} \rangle \left( 1 + \frac{1}{2} \frac{\sigma_{s_{\phi}}^2}{\langle s_{\phi} \rangle |\langle s_{\phi} \rangle|} \mathrm{Da}_{r_{\phi}} \right), \qquad (43) $$

where $\phi = w$ or $\phi = i$, as before. These approximations are derived in Appendix C for a simplified particle-growth model, and assuming small particle Damköhler numbers, small supersaturation Damköhler numbers, and narrow particle-size distributions. Equation (43) shows that the turbulence term $\langle r'_{\phi} s'_{\phi} \rangle$ is proportional to the particle Damköhler numbers. It can therefore be neglected when $\mathrm{Da}_{r_w}$ and $\mathrm{Da}_{r_i}$ are small, when the particles are large enough. Neglecting $\langle r'_{\phi} s'_{\phi} \rangle$ corresponds to the decoupling between the fluctuations of radii and supersaturation when $s'_{\phi}$ evolves much faster than $r'_{\phi}$ for large particles, so that $r'_{\phi}$ can respond only to the slow evolving $\langle s_{\phi} \rangle$, but not $s'_{\phi}$. We note that Eq. (43) fails near $\langle s_{\phi} \rangle \approx 0$, when the turbulence term $\langle r'_{\phi} s'_{\phi} \rangle$ becomes comparable to $\langle r_{\phi} \rangle \langle s_{\phi} \rangle$. In practice this affects most the stationary state where the system spends significant time at $s_i \approx 0$, i.e. after glaciation has happened.

Now consider the evolution of the average particle radii. In Appendix C we derive an approximate equation for $\langle r_{\phi} \rangle$ under the same assumptions used to derive Eq. (43)

$$ \frac{\mathrm{d}\langle r_{\phi} \rangle}{\mathrm{d}t} \approx A_{3,\phi} \frac{\langle s_{\phi} \rangle}{\langle r_{\phi} \rangle} \left( 1 - \frac{1}{2} \frac{\sigma_{s_{\phi}}^2}{\langle s_{\phi} \rangle |\langle s_{\phi} \rangle|} \mathrm{Da}_{r_{\phi}} \right). \qquad (44) $$

This result shows that small-scale turbulent fluctuations are negligible for small particle Damköhler numbers. The physical reason why the mean radius dynamics is unaffected by turbulent fluctuations of supersaturation is the same decoupling between the dynamics of $r'_{\phi}$ and $s'_{\phi}$ discussed in the previous paragraph. Equation (44) comes with the same caveat as above, it fails near $s_{\phi} \approx 0$. We stress that the statistical model correctly describes the mean radius even when $s_{\phi} \approx 0$.

In conclusion, Figs. 1 and 2, as well as Eqs. (43) and (44) show that turbulence has at best a weak effect upon the evolution of the mean particle size at small Damköhler numbers. This implies that small-scale turbulence neither
affects the WBF process nor the resulting glaciation time, defined as the time it takes for the ice-particle mass fraction IWC/(IWC+LWC+q_vρ_0) to reach 90% [7].

On the other hand, Figure 3(a) shows that small-scale turbulence has an effect upon the fluctuations $\sigma_{r_\phi}^2 = \langle r_\phi^{\prime 2} \rangle$ of the particle radii, as Chen et al. [6] concluded in their DNS study. Our model can explain this as follows. Equations (16), $dr_\phi/dt \sim A_{3,\phi}s_\phi/r_\phi$ imply that particles evaporate more rapidly the smaller they are. As a consequence, any small differences in radii of evaporating particles amplify in a deterministic fashion, explaining the rapid growth of the variance in Fig. 3. However, when the initial variance vanishes as for the CTGC, turbulent supersaturation fluctuations are required to initially widen the particle-size distribution, triggering the amplification. The widening of $\sigma_{r_\phi}$ happens for water droplets, but not for ice particles. Since the ice particles grow, the variance $\langle r_\phi^{\prime 2} \rangle$ shrinks instead. We note, finally, that the peaks of $\sigma_{r_w}/\langle r_w \rangle$ in Fig. 3 coincide with the peaks of the droplet Damköhler number $Da_{r_w}$ [Fig. 3(b)], confirming that larger particle Damköhler numbers imply stronger coupling between $s'_\phi$ and $r'_\phi$.

The above reasonings regarding the effect of small-scale turbulence on the average particle radii and their fluctuations have to be adjusted when discussing the core of the Pi chamber (Fig. 4), because we did not account for the particle injection and removal. In this case the particle-size distribution of larger droplets ($r_w > 3.5 \ \mu m$ as in Fig. 4) is dominated not by turbulence, but by injection of small particles and removal of larger ones, so small-scale turbulence has a weaker effect, compared with the CTGC. However, if we look at the PDF of $r_w$ for all droplet sizes [Fig. 5(a)], we see that turbulence has no influence only on cases with low ice injection rates. These cases correspond to larger particle sizes and smaller droplet Damköhler number $Da_{r_w}$. Once the ice injection rate increases and the cloud glaciates, only small or even unactivated droplets remain. Their $Da_{r_w}$ is larger and they are sensitive to supersaturation fluctuations, which widen droplet size distribution compared to the case without turbulence. The distributions of ice particle size [Fig. 5(b)] is unaffected by turbulence, because ice remains large and $Da_{r_i}$ is small.

At larger Damköhler numbers, small-scale turbulence could have a larger effect, in particular for spatially inhomo- geneous initial conditions [11]. Consider for example increasing the simulation-box size in the DNS to increase the size of the unresolved turbulent scales. This increases the Damköhler numbers (41), causing small-scale turbulence to be more important. In [6, 7] the linear size of the simulation domain was about 0.2 m, the Damköhler number $Da_{r_w}$ does not exceed 1.5 [Fig. 3(b)]. In order to get Damköhler numbers of order 15, where turbulence is expected to have a stronger effect [11], the linear size should be at least about 20 m (since $\tau_L \sim (L^2/\epsilon)^{1/3}$). But as discussed above, the statistical model may fail when the Damköhler numbers become too large. In this case more refined approximations for the condensation rates are needed [11]. It could also be of interest to formulate models aimed to describe the large-Da limit, as used in dense evaporating sprays [42].

Here we considered high number densities of ice particles, typical for deeper clouds. Polar stratus clouds tend to have lower ice-particle number densities, of the order $10^{-3} \ cm^{-3}$ [2]. In this case it is hard to reach the steady state with DNS (although it is possible for the statistical model). Figure 7 of Chen et al. [6] shows results for the initial growth of ice particles at lower ice-particle number density (the lowest value is $8 \times 10^{-3} \ cm^{-3}$). We performed our own DNSs for these cases, and find good agreement with the statistical-model results, and with the deterministic limit of the model (not shown).

Abade and Albuquerque [43] found that turbulence increases cloud glaciation times, and that ice particles and supercooled droplets experience different supersaturation fluctuations, in apparent contradiction with the DNS of [6], and with our statistical-model results. This discrepancy is explained by the approximations for condensation and deposition rates used in [43], where the influence of neighbouring particles (ice particles or water droplets) is not taken into account, hindering the WBF process.

Korolev and Milbrandt [44] speculated that ice particles and water droplets may be locally unmixed at small scales. In this case one expects the supersaturation distributions to be non-Gaussian [11], causing the statistical model to fail, in the form used here. Instead, one must rely on improved approximations of the supersaturation dynamics, such as mapping-closure approximations [9–11]. In this case we certainly expect the deterministic model to be inaccurate, in particular for the glaciation time. The question is by how much local unmixing delays glaciation, and how much turbulence contributes to the acceleration of it. This question also relates to a recent field campaign in which supercooled clouds were seeded with ice nuclei to initiate the WBF [45]. For these experiments, turbulence appears to be crucial to expose the newly nucleated ice crystals to water droplets [46].

## VI. CONCLUSIONS

We analysed the effect of small-scale turbulence on the glaciation process in mixed-phase clouds using a statistical model similar to those used to describe droplet evaporation in warm clouds. We found that the model describes DNS results for the glaciation process very well, for the parameters in [6] corresponding to a cloud-top generating cell and also for the parameters specified in [7] corresponding to the core of the Michigan Pi cloud chamber.
The statistical-model analysis shows that small-scale turbulence has an overall small effect on ice growth. Small-scale turbulence affects the evaporation of water droplets in the CTGC in that it initialises the growth of the droplet-size variance. In the core of the Pi chamber, for the parameters specified by Wang et al. [7], small-scale turbulence matters only close to the glaciation transition, where the water droplets are very small. When the corresponding time scale for droplet evaporation decreases so that it is of the same order of magnitude as the Lagrangian mixing time, then small-scale turbulence may have a stronger effect.

Our calculations show more generally that the effect of turbulence on glaciation is expected to be larger at larger Damköhler numbers, i.e., on larger spatial scales L. Since mixed-phase clouds exhibit a large range of vertical and horizontal length scales, covering shallow stratiform clouds of just hundred meters depth and a horizontal extent of several tens to hundreds of kilometers to deep convective clouds extending across the entire troposphere, the potential for small-scale turbulence to affect the glaciation process in these clouds needs to be acknowledged. This is especially true since most models used in the atmospheric sciences (from large-scale global circulation models to comparably high-resolution large-eddy simulation models) do not represent the effects of small-scale turbulence on cloud microphysical processes such as the WBF process on length scales smaller than at least few tens to a couple of hundred meters (i.e., their grid spacing).

A statistical model, such as the one presented here, may not only enable us to assess the effects of small-scale turbulence on the WBF process on larger length scales in a follow-up study, but could also be the basis of a parameterisation to represent the effects of small-scale turbulence in larger-scale models, which still struggle to represent mixed-phase clouds, and especially the coupling of cloud microphysics and turbulence [47].

*Acknowledgments*. BM, FH, and AP thank E. Bodenschatz for numerous discussions about the role of turbulence in atmospheric clouds. B. Mehlig and G. Sartnitsky were supported by Vetenskapsrådet (grant no. 2021-4452). F. Hoffmann was supported by the Emmy Noether program of the German Research Foundation (DFG) under Grant HO 6588/1-1. Gaetano Sardina was supported by Vetenskapsrådet (grant no. 2023-2026) and by ERC grant MixClouds 101126050 funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725.

[1] T Storelvmo and I Tan, "The Wegener-Bergeron-Findeisen process – its discovery and vital importance for weather and climate," Meteorologische Zeitschrift **24**, 455 (2015).

[2] H Morrison, G De Boer, G Feingold, J Harrington, MD Shupe, and K Sulia, "Resilience of persistent arctic mixed-phase clouds," Nature Geoscience **5**, 11–17 (2012).

[3] T Storelvmo, JE Kristjánsson, U Lohmann, T Iversen, A Kirkevåg, and Ø Seland, "Modeling of the Wegener-Bergeron-Findeisen process—implications for aerosol indirect effects," Environmental Research Letters **3**, 045001 (2008).

[4] AV Korolev, "Rates of phase transformations in mixed-phase clouds," Quarterly Journal of the Royal Meteorological Society **134**, 595–608 (2008).

[5] B Ervens, G Feingold, K Sulia, and J Harrington, "The impact of microphysical parameters, ice nucleation mode, and habit growth on the ice/liquid partitioning in mixed-phase arctic clouds," Journal of Geophysical Research: Atmospheres **116** (2011).

[6] S Chen, L Xue, S Tessendorf, K Ikeda, C Weeks, R Rasmussen, M Kunkel, D Blestrud, S Parkinson, M Meadows, and N Dawson, "Mixed-phase direct numerical simulation: ice growth in cloud-top generating cells," Atmospheric Chemistry and Physics **23**, 5217–5231 (2023).

[7] A Wang, S Krueger, and S Chen, "Pi chamber mixed-phase cloud simulation case," (2024), accessed: 2024-05-23.

[8] K Chang, J Bench, M Brege, W Cantrell, K Chandrakar, D Ciochetto, C Mazzoleni, LR Mazzoleni, D Niedermeier, and RA Shaw, "A laboratory facility to study gas–aerosol–cloud interactions in a turbulent environment: The π chamber," Bulletin of the American Meteorological Society **97**, 2343–2358 (2016).

[9] SB Pope, "PDF methods for turbulent reactive flows," Progress in Energy and Combustion Science **11**, 119–192 (1985).

[10] H Chen, S Chen, and RH Kraichnan, "Probability distribution of a stochastically advected scalar field," Phys. Rev. Lett. **63**, 2657 (1989).

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

[12] M Kulmala, Ü Rannik, EL Zapadinsky, and CF Clement, "The effect of saturation fluctuations on droplet growth," Journal of Aerosol Science **28**, 1395–1409 (1997).

[13] C Siewert, J Bec, and G Krstulovic, "Statistical steady state in turbulent droplet condensation," J. Fluid Mech. **810**, 254–280 (2017).

[14] J Fries, G Sardina, G Svensson, and B Mehlig, "Key parameters for droplet evaporation and mixing at the cloud edge," QJRMS **147**, 2160–2172 (2021).
[15] GC Abade, WW Grabowski, and H Pawlowska, "Broadening of cloud droplet spectra through eddy hopping: Turbulent entraining parcel simulations," JAS **75**, 3365-3379 (2017).

[16] G Sardina, S Poulain, L Brandt, and R Caballero, "Broadening of cloud droplet size spectra by stochastic condensation: Effects of mean updraft velocity and CCN activation," JAS **75**, 451-467 (2018).

[17] AV Korolev and IP Mazin, "Supersaturation of water vapor in clouds," Journal of the Atmospheric Sciences **60**, 2957-2974 (2003).

[18] M. Pinsky, A. Khain, and A. Korolev, "Theoretical analysis of liquid-ice interaction in the unsaturated environment with application to the problem of homogeneous mixing," Journal of the Atmospheric Sciences **75**, 1045-1062 (2018).

[19] 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).

[20] 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).

[21] J Bec, K Gustavsson, and B Mehlig, "Statistical models for the dynamics of heavy particles in turbulence," Annual Reviews of Fluid Mechanics **56**, 189-213 (2024).

[22] MD Petters and SM Kreidenweis, "A single parameter representation of hygroscopic growth and cloud condensation nucleus activity," Atmospheric Chemistry and Physics **7**, 1961-1971 (2007).

[23] Yefim L Kogan, "The simulation of a convective cloud in a 3-d model with explicit microphysics. part i: Model description and sensitivity experiments," Journal of the Atmospheric Sciences **48**, 1160-1189 (1991).

[24] CA Jeffery, "Inhomogeneous cloud evaporation, invariance, and Damköhler number," J. Geophys. Res. **112** (2007).

[25] W Mordy, "Computations of the growth by condensation of a population of cloud droplets," Tellus **11**, 16-44 (1959).

[26] PA Vaillancourt, MK Yau, P Bartello, and Wojciech W Grabowski, "Microscopic approach to cloud droplet growth by condensation. Part II: Turbulence, clustering, and condensational growth," JAS **59**, 3421-3435 (2002).

[27] R Paoli and K Shariff, "Turbulent condensation of droplets: Direct simulation and a stochastic model," JAS **66**, 723-740 (2009).

[28] I Saito, T Gotoh, and T Watanabe, "Broadening of cloud droplet size distributions by condensation in turbulence," Journal of the Meteorological Society of Japan. Ser. II **97**, 867-891 (2019).

[29] KK Chandrakar, W Cantrell, K Chang, D Ciochetto, D Niedermeier, M Ovchinnikov, RA Shaw, and F Yang, "Aerosol indirect effect from turbulence-induced broadening of cloud-droplet size distributions," PNAS **113**, 14243-14248 (2016).

[30] G Sarnitsky and S Heinz, "Nonparametric inference for diffusion processes in systems with smooth evolution," Physica A: Statistical Mechanics and its Applications **598**, 1273-1306 (2022).

[31] JR Ristorcelli, "Passive scalar mixing: Analytic study of time scale ratio, variance, and mix rate," Physics of Fluids **18** (2006).

[32] S Chen, "Replication Data for "Mixed-phase Direct Numerical Simulation: Ice Growth in Cloud-Top Generating Cells"," (2022).

[33] F Hoffmann, B Mayer, and G Feingold, "A parameterization of interstitial aerosol extinction and its application to marine cloud brightening," J. Atmos. Sci., 2849-2862 (2022).

[34] A Pumir, B Shraiman, and ED Siggia, "Exponential tails and random advection," Phys. Rev. Lett. **66**, 2984-2987 (1991).

[35] Z Warhaft et al., "Probability distribution of a passive scalar in grid-generated turbulence," Phys. Rev. Lett. **67**, 3503 (1991).

[36] JP Gollub, J Clarke, M Gharib, B Lane, and ON Mesquite, "Fluctuations and transport in a stirred fluid with a mean gradient," Phys. Rev. Lett. **67**, 3507 (1991).

[37] Jayesh and Z Warhaft, "Probability distribution, conditional dissipation, and transport of passive temperature fluctuations in grid-generated turbulence," Phys. Fluids A **4**, 2292-2307 (1992).

[38] Z Warhaft, "Passive scalars in turbulent flows," ARFM **32**, 203-240 (2000).

[39] B Kumar, F Janetzko, J Schumacher, and RA Shaw, "Extreme responses of a coupled scalar-particle system during turbulent mixing," NJP **14**, 115020 (2012).

[40] B Kumar, J Schumacher, and RA Shaw, "Lagrangian mixing dynamics at the cloudy-clear air interface," JAS **71**, 2564-2580 (2014).

[41] E Villermaux, "Mixing vs. stirring," Annu. Rev. Fluid Mech. **51**, 256-273 (2019).

[42] 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).

[43] GC Abade and DG Albuquerque, "Persistent mixed-phase states in adiabatic cloud parcels under idealised conditions," Quarterly Journal of the Royal Meteorological Society (2024).

[44] AV Korolev and J Milbrandt, "How are mixed-phase clouds mixed?" Geophysical Research Letters **49**, e2022GL099578 (2022).

[45] J Henneberger, F Ramelli, R Spirig, N Omanovic, AJ Miller, C Fuchs, H Zhang, J Bühl, M Hervo, ZA Kanji, et al., "Seeding of supercooled low stratus clouds with a uav to study microphysical ice processes: An introduction to the cloudlab project," Bull. Am. Meteorol. Soc. **104**, E1962-E1979 (2023).

[46] N Omanovic, S Ferrachat, C Fuchs, J Henneberger, AJ Miller, K Ohneiser, F Ramelli, P Seifert, R Spirig, H Zhang, et al., "Evaluating the Wegener-Bergeron-Findeisen process in icon in large-eddy mode with in situ observations from the cloudlab project," Atmos. Chem. Phys. **24**, 6825-6844 (2024).

[47] É Vignon, SP Alexander, PJ DeMott, G Sotiropoulou, F Gerber, TCJ Hill, R Marchand, A Nenes, and A Berne, "Challenging and improving the simulation of mid-level mixed-phase clouds over the high-latitude southern ocean," Journal of Geophysical Research: Atmospheres **126**, e2020JD033490 (2021).
[48] World Meteorological Organization, "Guide to instruments and methods of observation," (2021).

[49] MK Yau and RR Rogers, *A Short Course in Cloud Physics* (Butterworth Heinemann, 1989).

[50] D Lamb and J Verlinde, *Physics and chemistry of clouds* (Cambridge University Press, 2011).

[51] HR Pruppacher and JD Klett, "Microphysics of clouds and precipitation—kluwer academic publishers," Atmos. Oceanogr. Sci. Lib. **18**, 433–446 (1997).

[52] EL Andreas *et al.*, *Handbook of Physical Constants and Functions for Use in Atmospheric Boundary Layer Studies*, AD-a440 352 (Cold Regions Research and Engineering Laboratory (US, 2005).

[53] I Sölch and B Kärcher, "A large-eddy model for cirrus clouds with explicit aerosol and ice microphysics and Lagrangian ice particle tracking," Quarterly Journal of the Royal Meteorological Society **136**, 2074–2093 (2010).

[54] G Sarnitsky, *Inference methods for stochastic turbulence modeling with a closure perspective*, Ph.D. thesis, University of Wyoming (2021).

[55] WW Grabowski, M Andrejczuk, and L-P Wang, "Droplet growth in a bin warm-rain scheme with twomey ccn activation," Atmospheric Research **99**, 290–301 (2011).

[56] MK Chen, Sand Yau, P Bartello, and L Xue, "Bridging the condensation–collision size gap: a direct numerical simulation of continuous droplet growth in turbulent clouds," Atmospheric Chemistry and Physics **18**, 7251–7262 (2018).

[57] S Chen, L Xue, and MK Yau, "Impact of aerosols and turbulence on cloud droplet growth: an in-cloud seeding case study using a parcel–DNS (direct numerical simulation) approach," Atmospheric Chemistry and Physics **20**, 10111–10124 (2020).

[58] WJ Massman, "A review of the molecular diffusivities of H₂O, CO₂, CH₄, CO, O₃, SO₂, NH₃, N₂O, NO, and NO₂ in air, O₂ and N₂ near STP," Atmospheric environment **32**, 1111–1127 (1998).

[59] WJ Massman, "Molecular diffusivities of Hg vapor in air, O₂ and N₂ near STP and the kinematic viscosity and thermal diffusivity of air near STP," Atmospheric environment **33**, 453–457 (1999).

[60] J Wang, F Dalla Barba, A Roccon, G Sardina, A Soldati, and F Picano, "Modelling the direct virus exposure risk associated with respiratory events," Journal of the Royal Society Interface **19**, 20210819 (2022).

[61] SB Pope, *Turbulent Flows* (Cambridge University Press, 2000).
Appendix A: Model parameters

In this Appendix we summarise the details necessary to understand the definition and values of the model parameters used in the main text, for the CTGC [6] and for the Pi chamber [7].

*Saturation pressures.* The saturation pressures $p_{v,w}$ and $p_{v,i}$ are functions of only temperature. For the Pi chamber we take [48]

$$p_{v,w}(T) = 611.2 \text{ Pa} \exp\left(17.62 \frac{T - 273.15 \text{ K}}{T - 30.03 \text{ K}}\right), \qquad \text{(A1a)}$$

$$p_{v,i}(T) = 611.2 \text{ Pa} \exp\left(22.46 \frac{T - 273.15 \text{ K}}{T - 0.53 \text{ K}}\right). \qquad \text{(A1b)}$$

For CTGC we take (Chen, private communication)

$$p_{v,w}(T) = 2.53 \times 10^{11} \text{ Pa} \exp\left(-\frac{5.42 \times 10^3 \text{ K}}{T}\right), \qquad \text{(A2)}$$

$$p_{v,i}(T) = 3.41 \times 10^{12} \text{ Pa} \exp\left(-\frac{6.13 \times 10^3 \text{ K}}{T}\right). \qquad \text{(A3)}$$

*Parameters A<sub>3,φ</sub> and r<sub>A3,φ</sub>*. To derive expressions (19) and (18), we start from the well-known form of the particle growth equations [25, 49]:

$$\frac{dr_w^2}{dt} = 2\hat{A}_{3,w}(r_w) [s_w - s_{w,K}(r_w)], \qquad \text{(A4a)}$$

$$\frac{dr_i^2}{dt} = 2\hat{A}_{3,i}(r_i) a_3(r_i) s_i. \qquad \text{(A4b)}$$

To derive expressions (19) and (18) we need to show that $\hat{A}_{3,\phi}(r_\phi) = A_{3,\phi}a_3(r_\phi/r_{A_{3,\phi}})$. The radius-dependent function $\hat{A}_{3,\phi}$ in Eq. (A4) is given by

$$\hat{A}_{3,\phi} = \left[ \left( \frac{L_{\phi}(T_0)}{R_{v}T_0} - 1 \right) \frac{R_a}{c_p} \frac{\rho_{\phi} L_{\phi}(T_0)}{\kappa'_{T,\phi} p_0} + \frac{\rho_{\phi} R_{v} T_0}{\kappa'_{q_v,\phi} p_{v,\phi}(T_0)} \right]^{-1}, \qquad \text{(A5a)}$$

$$\kappa'_{T,\phi}(r_{\phi}) = \kappa_T \left( \frac{r_{\phi}}{r_{\phi} + \Delta_T} + \frac{\kappa_T}{r_{\phi} \alpha_{T,\phi}} \sqrt{\frac{2\pi}{R_a T_0}} \right)^{-1}, \qquad \text{(A5b)}$$

$$\kappa'_{q_v,\phi}(r_{\phi}) = \kappa_{q_v} \left( \frac{r_{\phi}}{r_{\phi} + \Delta_{q_v}} + \frac{\kappa_{q_v}}{r_{\phi} \alpha_{q_v,\phi}} \sqrt{\frac{2\pi}{R_v T_0}} \right)^{-1}. \qquad \text{(A5c)}$$

This form for $\hat{A}_{3,\phi}$ was used by [6] and [7], however we note that [50] provide a different formula that involves the Köhler correction function $s_{w,K}$. Here $\Delta_T$ and $\Delta_{q_v}$ are scales of the order of the mean free path in air, $\lambda_a$, parameterising kinetic corrections. The parameters $\alpha_{T,\phi}$ and $\alpha_{q_v,\phi}$ are thermal and condensation accommodation coefficients. We take their values from Refs. [51–53],

$$\begin{aligned} \Delta_T &= 2.16 \times 10^{-7} \text{ m}, & \Delta_{q_v} &= 0.87 \times 10^{-7} \text{ m}, \\ \alpha_{T,w} &= \alpha_{T,i} = 0.7, & \alpha_{q_v,w} &= 0.036, & \alpha_{q_v,i} &= 0.5. \end{aligned} \qquad \text{(A6)}$$

For these values, the dominant radius-dependent term in Eq. (A5a) is the one with $\kappa'_{q_v,\phi}$. Thus we the ignore kinetic corrections and the temperature accommodation correction, effectively setting $\Delta_T = \Delta_{q_v} = 0$ and $\kappa_{T,\phi} = \kappa_T$. In addition, we observe that $L_\phi/(R_v T_0) \gg 1$ and proceed to derive $\hat{A}_{3,\phi} = A_{3,\phi}a_3(r_\phi/r_{A_{3,\phi}})$:

$$\begin{aligned} \hat{A}_{3,\phi} &= \left[ \left( \frac{L_{\phi}(T_0)}{R_{v}T_0} - 1 \right) \frac{R_a}{c_p} \frac{\rho_{\phi} L_{\phi}(T_0)}{\kappa_{T,\phi} p_0} + \frac{\rho_{\phi} R_{v} T_0}{\kappa_{q_v,\phi} p_{v,\phi}(T_0)} \right]^{-1} \\ &= \left( \frac{R_a}{R_v} \frac{\rho_{\phi} L_{\phi}^2(T_0)}{c_p T_0 p_0} + \frac{\rho_{\phi} R_v T_0}{\kappa_{q_v} p_{v,\phi}(T_0)} + \frac{\rho_{\phi} \sqrt{2\pi R_v T_0}}{r_{\phi} \alpha_{q_v,\phi} p_{v,\phi}(T_0)} \right)^{-1} \qquad \text{(A7)} \\ &= \left( \frac{1}{A_{3,\phi}} + \frac{1}{A_{3,\phi}} \frac{r_{A_{3,\phi}}}{r_{\phi}} \right)^{-1} = A_3 a_3(r/r_{A_{3,\phi}}). \end{aligned}$$
The last row yields Eqs. (19) and (18).

*Supersaturation variance and correlation time.* Now we describe how we obtained $\sigma_{s_w}$ and $\tau_{s_w}^{(L)}$ for the CTGC [6] and for the Pi chamber [7].

The value of $\sigma_{s_w}$ is specified for the Pi chamber, but for the CTGC only the standard deviations $\sigma_{q_v}$ and $\sigma_T$ of the mixing ratio and temperature are provided. We reconstruct the value of $\sigma_{s_w}$ from the data provided in the supplementary material for [6]. Recall that $A_4(T) = p_{v,w}(T)/p_{v,i}(T)$, so expanding $A_4$ near $T = \langle T \rangle$ we express $s_i$ as

$$s_i = A_4(T)(s_w + 1) - 1 = \left[ A_4(\langle T \rangle) + T' \frac{dA_4}{dT} \bigg|_{T=\langle T \rangle} + \dots \right] (s_w + 1) - 1. \quad (A8)$$

After averaging this expression we can express the correlation $\langle s'_w T' \rangle$, Eq. (A9a). At the same time, we can compute the correlations $\langle s'_w T' \rangle$ and $\langle s'^2_w \rangle$ straight from Eq. (6a), taking $T_0 = \langle T \rangle$, $p_0 = \langle p \rangle$, $s_{w,0} = \langle s_w \rangle$. As before, we note that in our case the pressure term is negligible within the Oberbeck-Boussinesq approximation. Overall we get a system of three equations:

$$\langle s'_w T' \rangle = \left[ (\langle s_i \rangle + 1) - A_4(\langle T \rangle) (\langle s_w \rangle + 1) \right] / \left. \frac{dA_4}{dT} \right|_{T=\langle T \rangle}, \quad (A9a)$$

$$\langle s'_w T' \rangle = \frac{R_v}{R_a} \frac{\langle p \rangle}{p_{v,w}(\langle T \rangle)} \langle q'_v T' \rangle - (1 + \langle s_w \rangle) \frac{L_w(\langle T \rangle)}{R_v \langle T \rangle^2} \langle T'^2 \rangle, \quad (A9b)$$

$$\langle s_w'^2 \rangle = \left[ \frac{R_v}{R_a} \frac{\langle p \rangle}{p_{v,w}(\langle T \rangle)} \right]^2 \langle q_v'^2 \rangle + \left[ (1 + \langle s_w \rangle) \frac{L_w(\langle T \rangle)}{R_v \langle T \rangle^2} \right]^2 \langle T'^2 \rangle \\ - 2(1 + \langle s_w \rangle) \frac{R_v}{R_a} \frac{\langle p \rangle}{p_{v,w}(\langle T \rangle)} \frac{L_w(\langle T \rangle)}{R_v \langle T \rangle^2} \langle q_v' T' \rangle. \quad (A9c)$$

We solve this system to express $\sigma_{s_w}^2 = \langle s_w'^2 \rangle$ as a function of $\sigma_{q_v}^2 = \langle q_v'^2 \rangle$, $\sigma_T^2 = \langle T'^2 \rangle$, $\langle s_w \rangle$, $\langle s_i \rangle$, $\langle T \rangle$, $\langle p \rangle$; these datasets are provided in [32].

The calculation of $\tau_{s_w}$ as the parameter of the model (35) follows the procedure described in [54, Eq. 5.2 and Appendix F].

*Thermodynamic parameter values.* The thermodynamic parameters for the CTGC case are given in Table IV. They were taken from either directly from [6] and its Replication Data [32], or from the previous papers in the series [55–57]. For the reference value of supersaturation $s_{w,0}$ we take the mean of $s_w = 0$ (water vapour saturated w.r.t. water) and $s_w = 1/A_4 - 1$ (water vapour saturated w.r.t. ice):

$$s_{w,0} = \frac{1}{2} \left( \frac{1}{A_4} - 1 \right). \quad (A10)$$

This value of $s_{w,0}$ should be the a good general choice for describing WBF process in which supersaturation lies between these two points.

The parameter values for the Pi chamber [7] are summarised in Table V. For the diffusivities of temperature and water-vapour mixing ratio, $\kappa_T$ and $\kappa_{q_v}$, as well as the kinematic viscosity $\nu$ of air, we used [58, 59]:

$$\nu(p, T) = 1.327 \times 10^{-5} \frac{\mathrm{m}^2}{\mathrm{s}} \frac{101\,325\,\mathrm{Pa}}{\mathrm{p}} \left( \frac{T}{273.15\,\mathrm{K}} \right)^{1.81}, \quad (A11a)$$

$$\kappa_T(p, T) = 1.869 \times 10^{-5} \frac{\mathrm{m}^2}{\mathrm{s}} \frac{101\,325\,\mathrm{Pa}}{\mathrm{p}} \left( \frac{T}{273.15\,\mathrm{K}} \right)^{1.81}, \quad (A11b)$$

$$\kappa_{q_v}(p, T) = 2.178 \times 10^{-5} \frac{\mathrm{m}^2}{\mathrm{s}} \frac{101\,325\,\mathrm{Pa}}{\mathrm{p}} \left( \frac{T}{273.15\,\mathrm{K}} \right)^{1.81}. \quad (A11c)$$

The value for $s_{w,0}$ was calculated from (5) from the values of $q_{v,0}$ and $T_0$ specified by [60].

## Appendix B: Condensation rates in the statistical model

In this Appendix we derive the expressions for the conditional condensations rates $\langle C_w | s_w, t \rangle$ and $\langle C_i | s_w, t \rangle$, the mean condensation rates $\langle C_w \rangle(t)$ and $\langle C_i \rangle(t)$, and the statistical model (36). The notation $\langle \cdot | s_w, t \rangle$ is a shorthand for the usual notation for the conditional averages:

$$\langle \cdot | s_w, t \rangle = \langle \cdot | S_w(x, t) = s_w \rangle. \quad (B1)$$
This defines the ensemble average over all flow realisations for which the supersaturation field at point $x$ at time $t$ equals $s_w$. We use upper- and lower-case letters to distinguish the random variable $S(x, t)$ from its value $s_w$. All results here are valid for a statistically homogeneous system. For more details on the mathematical tools used here see Appendix H in Ref. [61].

*Conditional condensation rates.* To simplify the notation, we ignore here the subscript $w$ in $s_w$ and use just $s$ for supersaturation, especially since the derivation is valid for $s_i$ too. First we show that any conditional average for quantities of the form (24) can be computed as

$$ \left\langle \sum_{\alpha=1}^{N_{\phi}} G(x - x_{\alpha}) F(r_{\alpha}, S(x_{\alpha}, t)) \middle| S(x, t) = s \right\rangle = \frac{N_{\phi}}{V} \frac{f_{\phi}(s, t)}{f(s, t)} \langle F | s, t \rangle_{\phi}. \quad (B2) $$

Here $F$ is any function of a particle radius and supersaturation at the position of this particle, and $V$ is the volume of the simulation box. Next, $f(s, t)$ is the PDF of the Eulerian supersaturation field, it can be written as

$$ f(s, t) = \langle \delta(S(x, t) - s) \rangle, \quad (B3) $$

where $\delta$ is the Dirac delta function. The PDF $f_{\phi}$ of supersaturation at the particle positions (water droplets or ice particles, $\phi = w$ or $\phi = i$) is:

$$ f_{\phi}(s, t) = \frac{1}{N_{\phi}} \sum_{\alpha=1}^{N_{\phi}} \langle \delta(S(x_{\alpha}, t) - s) \rangle. \quad (B4) $$

Averages over the particle positions are denoted as $\langle \cdot \rangle_{\phi}$, $\phi = w$ or $\phi = i$:

$$ \langle F | s, t \rangle_{\phi} = \frac{1}{N_{\phi}} \sum_{\alpha=1}^{N_{\phi}} \langle F(r_{\alpha}, S(x_{\alpha}, t)) | S(x_{\alpha}, t) = s \rangle. \quad (B5) $$

To derive (B2), we use that the conditional average of a field $H(x, t)$ can be expressed as

$$ \langle H | S(x, t) = s \rangle = \frac{1}{f(s, t)} \langle H \delta(S(x, t) - s) \rangle. \quad (B6) $$

Second, we use that the system is statistically homogeneous. This allows us to take the spatial average of the l.h.s. of Eq. (B2), resulting in

$$ \left\langle \sum_{\alpha=1}^{N_{\phi}} G(x - x_{\alpha}) F(r_{\alpha}, S(x_{\alpha}, t)) \middle| S(x, t) = s \right\rangle = \sum_{\alpha=1}^{N_{\phi}} \frac{1}{f(s, t)} \langle G(x - x_{\alpha}) F(r_{\alpha}, S(x_{\alpha}, t)) \delta(S(x, t) - s) \rangle \quad (B7) \\ = \frac{1}{V} \int_V \frac{1}{f(s, t)} \sum_{\alpha=1}^{N_{\phi}} \langle G(x - x_{\alpha}) F(r_{\alpha}, S(x_{\alpha}, t)) \delta(S(x, t) - s) \rangle d\mathbf{x}. $$

Third, we require that the support of spatial kernels $G$ is smaller than the length scale at which $s(x, t)$ varies, so that we can treat $G$ as a delta function. Under these conditions we obtain

$$ \left\langle \sum_{\alpha=1}^{N_{\phi}} G(x - x_{\alpha}) F(r_{\alpha}, S(x_{\alpha}, t)) \middle| S(x, t) = s \right\rangle = \frac{1}{V} \frac{1}{f(s, t)} \sum_{\alpha=1}^{N_{\phi}} \langle F(r_{\alpha}, S(x_{\alpha}, t)) \delta(S(x_{\alpha}, t) - s) \rangle \quad (B8) \\ = \frac{1}{V} \frac{f_{\phi}(s, t)}{f(s, t)} \sum_{\alpha=1}^{N_{\phi}} \langle F(r_{\alpha}, S(x_{\alpha}, t)) | S(x_{\alpha}, t) = s \rangle = \frac{N_{\phi} f_{\phi}(s, t)}{V f(s, t)} \langle F | s, t \rangle_{\phi}. $$

Using this result and Eq. (23), we can rewrite the conditional condensation rates as

$$ \langle C_w | s_w, t \rangle = \frac{4}{3} \pi \frac{\rho_w}{\rho_0} \frac{N_w}{V} \frac{f_w(s_w, t)}{f(s_w, t)} \left\langle \frac{dr_w^3}{dt} \middle| s_w, t \right\rangle_w, \quad (B9a) $$

$$ \langle C_i | s_w, t \rangle = \frac{4}{3} \pi \frac{\rho_i}{\rho_0} \frac{N_i}{V} \frac{f_i(s_w, t)}{f(s_w, t)} \left\langle \frac{dr_i^3}{dt} \middle| s_w, t \right\rangle_i, \quad (B9b) $$
These formulae are consistent with Eq. (5) in [11], who considered a turbulent mixing problem with spatially inhomogeneous initial conditions and strong phase change, using mapping-closure approximations.

*Mean condensation rates.* Averaging expressions (B9) for $\langle C_w | s_w, t \rangle$ and $\langle C_i | s_w, t \rangle$ over $s_w$, we obtain Eqs. (33) for the mean condensation rates.

*Statistical model for $s'_w$ with condensation-rate fluctuations.* Now we derive the model (36) for the fluctuating supersaturation $s'_w$ that involves the condensation terms following the method of [30]. We start with the exact equation for $s'_w$, which for a statistically homogeneous system follows from (12a) as

$$ \frac{ds'_w}{dt} = \kappa \frac{\partial^2 s'_w}{\partial x_j \partial x_j} - A_{2,w} C'_w - A_{2,i} C'_i + f^{(s_w)'} \quad (B10) $$

If $s'_w$ has a fast oscillating component we can model it as a stochastic differential equation

$$ ds'_w = D^{(1)}(s'_w, t) dt + \sqrt{D^{(2)}(s'_w, t)} dW. \quad (B11) $$

This equation describes the behaviour of $s'_w$ on timescales larger than $\tau_{M,s_w}$, the Markov–Einstein timescale of $s_w$ and is not applicable for modelling the real behavior of $s'_w$ on smaller timescales. The Markov–Einstein timescale is of the order of the Taylor timescale for $s_w$, $\tau_{M,s_w} \sim \sqrt{\langle s_w'^2 \rangle / \langle (ds_w'/dt)^2 \rangle}$ and is of the order of the usual velocity Taylor timescale. In Itô's stochastic calculus that we use in Eq. (B11), future values of white noise are independent on the current value of $s'_w$ (the non-anticipation property). Since from (B11) we can express the white-noise term with

$$ \sqrt{D^{(2)}(s'_w, t)} dW = ds'_w - D^{(1)}(s'_w, t) dt \quad (B12) $$

we require it to obey the non-anticipation property:

$$ \left\langle \frac{ds'_w}{dt}(t + \Delta t) - D^{(1)}(s_w(t + \Delta t), t + \Delta t) \mid s_w, t \right\rangle = 0, \quad \text{for } \Delta t \ge \tau_{M,s_w}. \quad (B13) $$

Since $D^{(1)}$ evolves on timescales larger than $\tau_{M,s_w}$, we can approximate it from (B13) as

$$ D^{(1)}(s_w, t) = \left\langle \frac{ds'_w}{dt}(t + \tau_{M,s_w}) \mid s_w, t \right\rangle = \kappa \left\langle \frac{\partial^2 s'_w}{\partial x_j \partial x_j}(t + \tau_{M,s_w}) \mid s_w, t \right\rangle + \langle f^{(s_w)'}(t + \tau_{M,s_w}) \mid s_w \rangle \\ - A_{2,w} \langle C'_w(t + \tau_{M,s_w}) \mid s_w, t \rangle - A_{2,i} \langle C'_i(t + \tau_{M,s_w}) \mid s_w, t \rangle. \quad (B14) $$

Because the change of $s_w$ due to evaporation/condensation is slow compared to $\tau_{M,s_w}$ in our setup, we can substitute $\langle C'_w(t + \tau_{M,s_w}) \mid s_w, t \rangle$ with $\langle C'_w(t) \mid s_w, t \rangle$. For the diffusion term, we can use a usual Langevin mixing closure, while the forcing term will be accounted with $D^{(2)}$ since its only role is to keep $\sigma_{s_w}$ constant. Thus we turn (B11) into

$$ ds'_w = -A_{2,w} \langle C'_w | s'_w, t \rangle - A_{2,i} \langle C'_i | s'_w, t \rangle - \frac{1}{\tau_{s_w}^{(L)}} s'_w dt + \sqrt{D^{(2)}(s'_w, t)} dW(t). \quad (B15) $$

We could also relate $D^{(2)}(s_w, t)$ to the statistics of $ds'_w/dt$ as further described by [30]. However, for simplicity we just take $D^{(2)}$ to be independent of $s_w$, also ensuring that its value yields $\sigma_{s_w} = \text{const}$:

$$ D^{(2)} = \frac{2\sigma_{s_w}^2}{\tau_{s_w}^{(L)}} + 2A_{2,w}\langle C'_w s'_w \rangle + 2A_{2,i}\langle C'_i s'_w \rangle. \quad (B16) $$

Thus we arrive to Eq. (36).

## Appendix C: Effect of phase change and small-scale turbulence on particle-size distributions

In this Appendix we summarise details regarding the effect of small-scale turbulence on the mean supersaturation and the mean particle size, needed for the discussion in Section V. For simplicity we neglect radius-dependent corrections, replacing Eq. (16) by:

$$ \frac{dr_{\phi}^2}{dt} = \begin{cases} 2A_{3,\phi}s_{\phi} & \text{if } r_{\phi} > 0, \\ 0 & \text{if } r_{\phi} = 0. \end{cases} \quad (C1) $$
Further, we consider the state when none of the particle considered has evaporated completely. Also, the derivations below rely on expanding $1/r_φ$ around $1/\langle r_φ \rangle$ assuming small $r'_φ/\langle r_φ \rangle$,

$$ \frac{1}{r_{\phi}} = \frac{1}{\langle r_{\phi} \rangle} - \frac{r'_{\phi}}{\langle r_{\phi} \rangle^2} + \dots \qquad (C2) $$

and thus are valid for relatively sharp particle size distributions.

First we show that the correlation $\langle r'_φ s'_φ \rangle$ of the fluctuating quantities $r'_φ$ and $s'_φ$ is of the order of $\text{D}\alpha_{r_φ}$,

$$ \langle r'_{\phi} s'_{\phi} \rangle = \frac{1}{2} \text{D}\alpha_{r_{\phi}} \frac{\sigma_{s_{\phi}}^2}{\langle s_{\phi} \rangle |\langle s_{\phi} \rangle|} \langle r_{\phi} \rangle \langle s_{\phi} \rangle, \qquad (C3) $$

We start by deriving the evolution equation for $\langle r'_φ s'_φ \rangle$. From Eqs. (C1), (12b), (35) and (C2) we find

$$ \begin{aligned} \frac{d\langle r'_{\phi} s'_{\phi} \rangle}{dt} &= \frac{d\langle r_{\phi} s'_{\phi} \rangle}{dt} = \left\langle \frac{d(r_{\phi} s'_{\phi})}{dt} \right\rangle = \left\langle r_{\phi} \frac{ds'_{\phi}}{dt} \right\rangle + \left\langle s'_{\phi} \frac{dr_{\phi}}{dt} \right\rangle \\ &= -\frac{1}{\tau_{s_w}^{(L)}} \langle r'_{\phi} s'_{\phi} \rangle + A_{3,\phi} \left\langle \frac{s_{\phi} s'_{\phi}}{r_{\phi}} \right\rangle = -\frac{1}{\tau_{s_w}^{(L)}} \langle r'_{\phi} s'_{\phi} \rangle + A_{3,\phi} \frac{\sigma_{s_{\phi}}^2}{\langle r_{\phi} \rangle} + \dots \end{aligned} \qquad (C4) $$

For small Damköhler numbers, the evolution of the mean quantities, including $\langle r'_φ s'_φ \rangle$, happens on timescales much larger than the turbulent timescale $\tau_{s_w}^L$. Hence we can neglect the term $d\langle r'_{\phi} s'_{\phi} \rangle/dt$ compared to the term $\langle r'_{\phi} s'_{\phi} \rangle/\tau_{s_w}^{(L)}$. This results in

$$ \langle r'_{\phi} s'_{\phi} \rangle = \tau_{s_w}^{(L)} A_{3,\phi} \frac{\sigma_{s_{\phi}}^2}{\langle r_{\phi} \rangle}. \qquad (C5) $$

Using the definitions (42) and (41) of $\tau_{r_φ}$ and $\text{D}\alpha_{r_φ}$ we arrive at Eq. (C3).

Next, we derive (43). The mean condensation rates are

$$ \langle C_{\phi} \rangle = 4\pi \frac{\rho_{\phi}}{\rho_0} A_{3,\phi} n_{\phi} \langle r_{\phi} s_{\phi} \rangle = 4\pi \frac{\rho_{\phi}}{\rho_0} A_{3,\phi} n_{\phi} [\langle r_{\phi} \rangle \langle s_{\phi} \rangle + \langle r'_{\phi} s'_{\phi} \rangle]. \qquad (C6) $$

Using (C3) we obtain (43). The derivation of Eqs. (44) proceeds by averaging Eq. (16), using Eq. (C2) and Eq. (C3):

$$ \begin{aligned} \frac{d\langle r_{\phi} \rangle}{dt} &= A_{3,\phi} \left\langle \frac{s_{\phi}}{r_{\phi}} \right\rangle = A_{3,\phi} \left\langle \frac{s_{\phi}}{\langle r_{\phi} \rangle} - \frac{s_{\phi} r'_{\phi}}{\langle r_{\phi} \rangle} \right\rangle \qquad (C7) \\ &= A_{3,\phi} \frac{\langle s_{\phi} \rangle}{\langle r_{\phi} \rangle} \left( 1 - \frac{\langle r'_{\phi} s_{\phi} \rangle}{\langle r_{\phi} \rangle \langle s_{\phi} \rangle} \right) = A_{3,\phi} \frac{\langle s_{\phi} \rangle}{\langle r_{\phi} \rangle} \left( 1 - \frac{1}{2} \frac{\sigma_{s_{\phi}}^2}{\langle s_{\phi} \rangle |\langle s_{\phi} \rangle|} \text{D}\alpha_{r_{\phi}} \right). \end{aligned} $$
<table><caption>TABLE IV. Model parameters for Pi chamber [7].</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R<sub>a</sub></td><td>287.05 J/(kg K)</td></tr><tr><td>R<sub>v</sub></td><td>461.52 J/(kg K),</td></tr><tr><td>c<sub>p</sub></td><td>1005 J/(kg K)</td></tr><tr><td>T<sub>0</sub></td><td>265.63 K</td></tr><tr><td>p<sub>0</sub></td><td>1 &times; 10<sup>5</sup> Pa</td></tr><tr><td>&rho;<sub>0</sub></td><td>1.311 kg/m<sup>3</sup></td></tr><tr><td>V</td><td>8 &times; 10<sup>-3</sup> m<sup>3</sup></td></tr><tr><td>s<sub>w,0</sub></td><td>5.253 &times; 10<sup>-2</sup></td></tr><tr><td>q<sub>v,0</sub></td><td>2.28 &times; 10<sup>-3</sup> kg/kg</td></tr><tr><td>&nu;</td><td>1.278 &times; 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>&chi;<sub>T</sub></td><td>1.800 &times; 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>&chi;<sub>q<sub>v</sub></sub></td><td>2.098 &times; 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>&chi;</td><td>1.944 &times; 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>&rho;<sub>0</sub></td><td>1.311 kg/m<sup>3</sup></td></tr><tr><td>&rho;<sub>w</sub></td><td>1000 kg/m<sup>3</sup></td></tr><tr><td>&rho;<sub>i</sub></td><td>917 kg/m<sup>3</sup></td></tr><tr><td>&alpha;<sub>q<sub>v</sub>,w</sub></td><td>0.036</td></tr><tr><td>&alpha;<sub>q<sub>v</sub>,i</sub></td><td>0.5</td></tr><tr><td>r<sub>A<sub>3,w</sub></sub></td><td>2.805 &times; 10<sup>-6</sup> m</td></tr><tr><td>r<sub>A<sub>3,i</sub></sub></td><td>0.1906 &times; 10<sup>-6</sup> m</td></tr><tr><td>r<sub>dry</sub></td><td>0.0625 &times; 10<sup>-6</sup> m</td></tr><tr><td>r<sub>w,initial</sub></td><td>0.0625 &times; 10<sup>-6</sup> m</td></tr><tr><td>r<sub>i,initial</sub></td><td>2 &times; 10<sup>-6</sup> m</td></tr><tr><td>&kappa;</td><td>1.12</td></tr><tr><td>A<sub>2,w</sub></td><td>664.8</td></tr><tr><td>A<sub>2,i</sub></td><td>691.0</td></tr><tr><td>A<sub>3,w</sub></td><td>40.08 &times; 10<sup>-12</sup> m<sup>2</sup>/s</td></tr><tr><td>A<sub>3,i</sub></td><td>38.27 &times; 10<sup>-12</sup> m<sup>2</sup>/s</td></tr><tr><td>A<sub>4</sub></td><td>1.078</td></tr><tr><td>k<sub>&infin;,w</sub></td><td>1.233 &times; 10<sup>8</sup> m<sup>-1</sup>s<sup>-1</sup></td></tr><tr><td>k<sub>&infin;,i</sub></td><td>1.131 &times; 10<sup>8</sup> m<sup>-1</sup>s<sup>-1</sup></td></tr><tr><td>H</td><td>0.2 m</td></tr><tr><td>I<sub>w</sub></td><td>10 &times; <sup>10<sup>6</sup></sup>&frasl;<sub>60</sub> m<sup>-3</sup>s<sup>-1</sup></td></tr><tr><td>&sigma;<sub>T</sub></td><td>0.72 K</td></tr><tr><td>&sigma;<sub>q<sub>v</sub></sub></td><td>0.165 &times; 10<sup>-3</sup></td></tr><tr><td>&sigma;<sub>s<sub>w</sub></sub></td><td>2.047 &times; 10<sup>-2</sup></td></tr><tr><td>&tau;<sub>s<sub>w</sub></sub><sup>(L)</sup></td><td>0.755 s</td></tr><tr><td>&tau;<sub>s<sub>w</sub>,force</sub></td><td>60 s</td></tr><tr><td>s<sub>w,force</sub></td><td>5.253 &times; 10<sup>-2</sup></td></tr></tbody></table>
<table><caption>TABLE V. Model parameters for CTGC [6].</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R<sub>a</sub></td><td>287 J/(kg K)</td></tr><tr><td>R<sub>v</sub></td><td>467 J/(kg K),</td></tr><tr><td>c<sub>p</sub></td><td>1005 J/(kg K)</td></tr><tr><td>T<sub>0</sub></td><td>259.53 K</td></tr><tr><td>p<sub>0</sub></td><td>57 160 Pa</td></tr><tr><td>ρ<sub>0</sub></td><td>0.7674 kg/m<sup>3</sup></td></tr><tr><td>V</td><td>8 × 10<sup>-3</sup> m<sup>3</sup></td></tr><tr><td>s<sub>w,0</sub></td><td>-6.298 × 10<sup>-2</sup></td></tr><tr><td>q<sub>v,0</sub></td><td>2.171 × 10<sup>-3</sup> kg/kg</td></tr><tr><td>ν</td><td>1.6 × 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>κ<sub>T</sub></td><td>2.22 × 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>κ<sub>q<sub>v</sub></sub></td><td>2.55 × 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>κ</td><td>2.379 × 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>ρ<sub>w</sub></td><td>1000 kg/m<sup>3</sup></td></tr><tr><td>ρ<sub>i</sub></td><td>917 kg/m<sup>3</sup></td></tr><tr><td>α<sub>q<sub>v,w</sub></sub></td><td>0.036</td></tr><tr><td>α<sub>q<sub>v,i</sub></sub></td><td>0.036</td></tr><tr><td>r<sub>A<sub>3,w</sub></sub></td><td>3.313 × 10<sup>-6</sup> m</td></tr><tr><td>r<sub>A<sub>3,i</sub></sub></td><td>3.182 × 10<sup>-6</sup> m</td></tr><tr><td>r<sub>dry</sub></td><td>1 × 10<sup>-6</sup> m</td></tr><tr><td>κ</td><td>0.3</td></tr><tr><td>A<sub>2,w</sub></td><td>621.6</td></tr><tr><td>A<sub>2,i</sub></td><td>646.5</td></tr><tr><td>A<sub>3,w</sub></td><td>2.945 × 10<sup>-11</sup> m<sup>2</sup>/s</td></tr><tr><td>A<sub>3,i</sub></td><td>2.696 × 10<sup>-11</sup> m<sup>2</sup>/s</td></tr><tr><td>A<sub>4</sub></td><td>1.144</td></tr><tr><td>σ<sub>T</sub></td><td>0.143 K</td></tr><tr><td>σ<sub>q<sub>v</sub></sub></td><td>4.5 × 10<sup>-5</sup> kg/kg</td></tr><tr><td>τ<sub>s<sub>w</sub></sub><sup>(L)</sup></td><td>2.04 s</td></tr></tbody></table>