# Steady-state supersaturation distributions for clouds under turbulent
forcing

Manuel Santos Gutiérrez Corresponding author. Email: manuel.santos-
<EMAIL_ADDRESS>Department of Earth and Planetary Sciences, Weizmann
Institute of Science, Rehovot, Israel  Kalli Furtado Email:
<EMAIL_ADDRESS>Centre for Climate Research Singapore, Meteorological
Service Singapore, Singapore

###### Abstract

The supersaturation equation for a vertically moving adiabatic cloud parcel is
analysed. The effects of turbulent updrafts are incorporated in the shape of a
stochastic Lagrangian model, with spatial and time correlations expressed in
terms of turbulent kinetic energy. Using the Fokker-Planck equation, the
steady-state probability distributions of supersaturation are analytically
computed for a number of approximations involving the timescale separation
between updraft fluctuations and phase-relaxation, multiplicative noise
simplification and droplet or ice particle size fluctuations. While the
analytical results are presented in general for single-phase clouds, the
calculated distributions are used to compute mixed-phase cloud properties—
mixed fraction and mean liquid water content in an initially icy cloud— and
are argued to be useful for generalising and constructing new parametrisation
schemes.

###### Contents

  1. 1 Introduction
  2. 2 The supersaturation equation
    1. 2.1 The Squires equation
    2. 2.2 Stochastic model for updrafts
  3. 3 Quasi-steady model statistics
    1. 3.1 Fast decorrelation timescale: $\tau_{d}\ll 1$
      1. 3.1.1 Linear approximation
    2. 3.2 Slow decorrelation timescale: $\tau_{d}\gg 1$
  4. 4 Diffusional growth and size fluctuations
    1. 4.1 Fluctuations in droplet radius
  5. 5 Parametrisation of mixed-phase clouds
  6. 6 Discussion
  7. A List of some used notations and symbols
  8. B Stationary supersaturation distribution
    1. B.1 The general case
  9. C The square-root process

##  1 Introduction

Water vapour supersaturation is a key thermodynamic parameter in the formation
and development of warm and cold clouds. The activation of cloud condensation
nuclei (CCN), the spectrum of droplet sizes and, therefore, precipitation
heavily depend on the supersaturation budget of an evolving cloud. By
extension, the study of supersaturation is fundamental to determine the
radiative properties of cloud fields which constitute a major climate feedback
[10], yet a major source of uncertainty in climate projections [26].

Due to its microphysical character, supersaturation cannot be explicitly
resolved by global circulation models (GCMs) and, hence, it has to be
prescribed using adequate parametrisation techniques; see the seminal work in
[1]. One way of doing this is by predicting the probability of encountering
the cloud parcel at a certain relative humidity level in terms of the
background dynamical and large-scale information. Taking the GCM perspective,
the idea is to infer the subgrid cloud properties out of the prognosed
variables: vertical velocity, turbulent kinetic energy, temperature, humidity
and pressure, which altogether depict the relative humidity configuration of
an homogeneous cloud parcel. However, localised temperature gradients or
turbulence at the smaller scales can create inhomogeneous fluctuations which
alter the spread of supersaturation values, hence changing the microphysical
properties of the subgrid cloud parcel. This is particularly relevant in
mixed-phase clouds, where the turbulent dynamics influence the activation of
long-lived supercooled water in icy clouds [12, 9, 7].

Different updraft profiles— ranging from steady to turbulent— yield different
supersaturation distributions which determine key features like mean liquid
water content, mixed-phase cloud fraction or the activation fraction of CCN.
It is crucial, then, to determine not only the mean value of supersaturation
at given location, but also its variance or even higher moments. To this end,
the equations for supersaturation or condensational growth are coupled with a
suitable stochastic forcing law that captures the effects of turbulence [2,
21]. Hence, the classical deterministic models become stochastic differential
equations (SDEs), which have been widely used in the cloud physics community
[21, 10]. In fact, this stochastic physics framework has been used in the
elaboration of analytically tractable parametrisation schemes for mixed-phase
clouds [7], the study of droplet growth by condensation [2, 20] and the
determination of steady-state warm cloud properties [23]. Departing from the
work of [7], we extend their analytical predictions of supersaturation
distribution to a wider range of contexts, involving the equations of
turbulent updrafts, supersaturation and diffusional growth.

This paper is structured as follows. In Section 2, the Squires equation for
the evolution of supersaturation is revisited from first principles and,
secondly, a stochastic equation for turbulent updrafts is presented, in the
lines of the theory of stochastic Lagrangian turbulent models; see, e.g.,
[18]. In Section 3, the quasi-steady equation— which assumes a constant mean
droplet/particle radius— is analysed in a number of approximations providing
formulas for the probability distribution of supersaturation. In Section 4,
fluctuations in droplet size are allowed and their net effects on
supersaturation distribution are investigated in analytical terms. A total of
five different supersaturation distributions are obtained on analytical
grounds. In Section 5, the relevance of the five obtained probability density
functions is discussed and compared in the context of mixed-phase clouds.
Finally, a discussion over the results is done in Section 6. To supplement the
information in the main text, appendices are included to discuss some
technical topics related to the analysis of stochastic differential equations.

##  2 The supersaturation equation

###  2.1 The Squires equation

We consider a vertically moving cloud parcel containing a monodisperse family
of liquid water droplets or ice particles that are spatially uniformly
distributed. Furthermore, it is assumed that droplet/particle number remains
constant in time. The evolution of supersaturation in the cloud is,
essentially, determined by the sources and sinks of relative humidity due to
the adiabatic cooling of the parcel in ascent, and condensation of water
vapour onto the existing droplets’ or particles’ surface [19]. However, the
exact relation between the rate of change of supersaturation and its sinks and
sources is obtained by taking its derivative with respect to to time. We
recall that supersaturation, for either ice or water— is defined as:

$\mathrm{S}=\frac{e-E}{E},$ (1)

where $e$ is the water vapour pressure and $E$ the same although at saturation
over a flat surface of liquid water or ice. We shall not specify now whether
we are dealing with water or ice supersaturation because the stochastic
analysis will be done independently. It is noted, however, that the
calculations immediately below can be done for single or mixed-phase clouds;
see [13].

Taking the time derivative of $\mathrm{S}$ and employing the mass
conservation, temperature, Clausius-Clayperon and the quasi hydrostatic
approximation, P. Squires derived in 1952 an equation to describe the
evolution of the supersaturation budget [24]— see also the Appendix in [13]:

$\frac{1}{1+\mathrm{S}}\mathrm{d}\mathrm{S}=au-b\frac{\mathrm{d}q}{\mathrm{d}t},$
(2)

where $u$ and $q$ are the vertical velocity and liquid-water or ice mixing
ratio, respectively. The evolution of $\mathrm{S}$ obeys an equation with a
nonlinear term stemming from the time and temperature dependence of the
equilibrium water vapour pressure $E$ in Eq. (1). The equation (2) reveals
that only to leading order in small values of $\mathrm{S}$ do we obtain a
linear dependence on vertical velocity $u$ and water vapour
condensation/deposition, $\mathrm{d}q/\mathrm{d}t$. Such is the case of warm
clouds, where supersaturation levels do not typically exceed $2\%$ [15].
However, mixed-phase conditions arise precisely when ice supersaturation
fluctuates fully icy clouds spreading beyond small values [7], making the
nonlinearity in Eq. (2) more relevant. This will be discussed in the next
section.

In order to find a closed model for supersaturation, it is necessary to
include the equation describing vapour condensation/deposition,
$\mathrm{d}q/\mathrm{d}t$. At a fixed time, the mixing ratio depends on the
concentration, size and density of the particles in the following fashion
[19]:

$q=\frac{4\pi
N}{3\rho_{a}}\int_{0}^{1}\int_{-\infty}^{\infty}\int_{0}^{\infty}f(r,\rho,c)\rho
r^{3}\mathrm{d}r\mathrm{d}\rho\mathrm{d},$ (3)

where $f$ is the particle size distribution and the different variables are
understood for either liquid droplets or ice particles. It is worth noting
that the first integral ranging from zero to unity refers to the possible
distribution of capacitances, which reflects the efficiency with which
droplets/ice particles collect water vapour [19]. While the capacitance of
droplets is approximately $1$— since all of them are almost spheres—, ice
particles possess a wider range of possible shapes that are displayed under
certain temperature and humidity conditions, affecting the collection of water
vapour. Taking the time derivative of Eq. (3):

$\frac{\mathrm{d}q}{\mathrm{d}t}=\frac{4\pi
N}{\rho_{a}}\int_{0}^{1}\int_{-\infty}^{\infty}\int_{0}^{\infty}f(r,\rho,c)\rho
r^{2}\frac{\mathrm{d}r}{\mathrm{d}t}\mathrm{d}r\mathrm{d}\rho\mathrm{d}c.$ (4)

We now employ the diffusional growth law for droplets/particles, which states
that water vapour condensates/deposits proportionally to the available
supersaturation and inversely proportionally to the radius of the particle
[19]:

$\frac{\mathrm{d}r}{\mathrm{d}t}=\frac{cA_{\circ}\mathrm{S}}{r}.$ (5)

Since we are considering a monodisperse particle size and shape, the previous
equation does not only model the diffusional growth of a single particle but
of the whole population. Inserting Eq. (5) into Eq. (4) we find:

$\frac{\mathrm{d}q}{\mathrm{d}t}=\mathrm{S}\frac{4\pi c\rho
A_{\circ}Nr}{\rho_{a}}.$ (6)

Finally, to obtain the full evolution equation for $\mathrm{S}$, the
expression for water vapour absorption in Eq. (6) is substituted into Eq. (2).
The updraft term is not yet specified, although we anticipate that an SDE for
updrafts will be coupled to the supersaturation equation in Section 2.2.

It is clear, then, that a closed equation for the evolution of supersaturation
with respect to ice can be obtained by simply integrating Eq. (6) and plugging
into Eq. (2) to give:

$\frac{1}{1+\mathrm{S}}\mathrm{d}\mathrm{S}=-b\mathrm{S}(t)\sqrt{r(0)^{2}+2cA_{\circ}\int_{0}^{t}\mathrm{S}(s)\mathrm{d}s}\mathrm{d}t+au.$
(7)

The solution of this equation has been shown to converge to a steady state
value, which is understood as an equilibrium supersaturation, whereby the
changes in relative humidity are balanced by the absorption/release of water
vapour by the droplets or ice particles [13]. Numerical solutions of this
equation can be achieved, although its analytical treatment entails expansions
and approximations [5]. Equation (7) is a nonlinear integro-differential
equation, since it possesses a square-root of a memory or integral term which
accounts for all past supersaturation configurations. The memory term is, in
principle, only valid for all the time that a droplets grow without
evaporating, sedimenting or precipitating. If any of these processes occurs,
memory is lost and the Eq. (7) would have to be reinitialised. To account for
this issue, in Section 4 we shall develop a theoretical framework to
understand the net effects of the fluctuations in droplet or particle radius
on the evolution of supersaturation.

###  2.2 Stochastic model for updrafts

The aim of this section is to construct a stochastic model for updraft
fluctuations. For this, a cloudy parcel is supposed to be embedded in a
turbulent environment with a prescribed surrounding supersaturation
$\mathrm{S}_{E}$, which is allowed to be negative in case of dry, subsaturated
air. The vertical velocity $u$ in such cloud parcel is assumed to be
decomposed into its mean $\bar{u}$ and fluctuating part $u^{\prime}$ so that
$u=\bar{u}+u^{\prime}$. Given the small spatial scales considered herein,
turbulence is taken to be isotropic so that the variance of $u^{\prime}$,
$\sigma_{u}$, is related to the turbulent kinetic energy $\mathrm{TKE}$ of the
background flow in the following way:

$\mathrm{TKE}=\frac{1}{2}\left(\sigma^{2}_{x}+\sigma^{2}_{y}+\sigma^{2}_{u}\right)=\frac{3}{2}\sigma^{2}_{u}.$
(8)

The updraft variance, $\sigma_{u}^{2}$, together with the eddy dissipation
rate $\varepsilon$ provide an estimate of the exponential rate, $1/\tau_{d}$,
at which the fluctuations $u^{\prime}$ decorrelate in time. Such number is
called the Lagrangian decorrelation timescale [18]. Roughly speaking
$\tau_{d}$ indicates the amount of time needed for a turbulent flow to forget
its original configuration:

$\displaystyle\tau_{d}=\frac{2\sigma^{2}_{u}}{\epsilon C_{0}}.$ (9)

Hence, if we assume that vertical motion is homogeneous, random, with
stationary mean $\bar{u}$, variance $\sigma^{2}_{w}$ and decorrelation time
$\tau_{d}$, the simplest model for the vertical velocity is the following red
noise equation [18, §3.5]:

$\mathrm{d}u(t)=-\frac{1}{\tau_{d}}\left(u(t)-\bar{u}\right)\mathrm{d}t+\sqrt{\frac{2}{\tau_{d}}}\sigma_{u}\mathrm{d}W_{t},$
(10)

where $W_{t}$ denotes a standard Wiener process, which accounts for the
acceleration increments over $\mathrm{d}t$ time-units, that come from random
pressure fluctuations which decorrelate instantly on time. Equation (10) has
the structure of an Ornstein-Uhlenbeck (OU) process, for which there exists a
vast collection of analytical results [25, 14]. In particular, this one-
dimensional Gaussian process satisfies the following mean, variance and
correlation properties:

$\displaystyle\mathbb{E}\left[u(t)\right]$
$\displaystyle=e^{-t/\tau_{d}}u(0)+(1-e^{-t/\tau_{d}})\bar{u}$ (11a)
$\displaystyle\mathrm{Var}\left(u(t)\right)$
$\displaystyle=\sigma_{u}^{2}\left(1-e^{-2t/\tau_{d}}\right)$ (11b)
$\displaystyle\mathbb{E}\left[u(t)u(0)\right]$
$\displaystyle=\sigma_{u}^{2}e^{-t/\tau_{d}}.$ (11c)

where $u(0)$ is the initial value of the vertical velocity. Thus, as $t$ tends
to infinity, the solution of Eq. (10) will distribute according to a Gaussian
function with mean $\bar{u}$ and variance $\sigma^{2}_{u}$.

As a result of considering a cloud in contact with the exterior, turbulent
motion will also be in charge of mixing mass at the cloud edges; see [6] for a
recent discussion of cloud edges and transition zones in cumulus clouds in
terms of adiabaticity of their components. This process is assumed to be
acting at a constant rate and will drive the supersaturation at the cloud
towards an equilibrium value $\mathrm{S}_{E}$. In other words, turbulent
mixing rates determine the characteristic time $\tau_{mix}$ to homogenise a
cloud volume with is surrounding reservoir [10]:

$\tau_{mix}=\left(\frac{L^{2}}{\epsilon}\right)^{1/3},$ (12)

where $L$ is the characteristic length of the turbulent zone. Therefore, a
more general form for the supersaturation equation is studied so that mixing
at the cloud edges is also taken into account:

$\frac{1}{1+\mathrm{S}}\frac{\mathrm{d}\mathrm{S}}{\mathrm{d}t}=au-b\mathrm{S}r-\frac{1}{\tau_{mix}}\left(\mathrm{S}-\mathrm{S}_{E}\right).$
(13)

Finally, it is enough to couple this equation to Eq. (10) and the diffusional
growth equation (5) to obtain a closed set of equations for the time evolution
of supersaturation.

##  3 Quasi-steady model statistics

The diffusional growth rate is inversely proportional to the radius of the
droplet or particle in question. Therefore, small droplets or particles grow
faster compared to larger ones. In this sense, if Eq. (5) is initialised with
a large radius $r^{2}(0)$, it is expected that it will remain almost constant,
at least for the interval where the following inequality is satisfied; see
also [13, 5]:

$r^{2}(0)\gg\left|2cA_{\circ}\int_{0}^{t}\mathrm{S}(s)\mathrm{d}s\right|.$
(14)

Hence, assuming that changes in the size of the cloud droplets can be
neglected, we can take $r$ in Eq. (4) to be constant, equal to $\bar{r}$. This
is called the quasi-steady approximation [13]. The resulting equation loses
the time-dependence of the variable $r$ and reads as:

$\displaystyle\frac{\mathrm{d}\mathrm{S}}{\mathrm{d}t}$
$\displaystyle=-(B+C)\left(\mathrm{S}-\frac{C\mathrm{S}_{E}}{B+C}\right)(1+\mathrm{S})+a(1+\mathrm{S})u$
(15a) $\displaystyle\mathrm{d}u$
$\displaystyle=-\frac{1}{\tau_{d}}\left(u-\bar{u}\right)\mathrm{d}t+\sqrt{\frac{2}{\tau_{d}}}\sigma_{u}\mathrm{d}W_{t},$
(15b)

where two constants have been introduced:

$\displaystyle B$ $\displaystyle=bB_{0}N\bar{r};$ (16a) $\displaystyle C$
$\displaystyle=\left(\frac{\epsilon}{L^{2}}\right)^{1/3}.$ (16b)

Note that the radius $r$ is no longer a variable and that is absorbed into the
constant $B$. The target of this section is to derive analytically the
stationary statistics of model (15) in a variety of approximations that are
detailed in the subsections bellow.

###  3.1 Fast decorrelation timescale: $\tau_{d}\ll 1$

There are two different timescales involved in Eq. (15). One is set by the
constant $B+C$, as the characteristic time for the absorption of water vapour,
and the other is the Lagrangian decorrelation timescale $\tau_{d}$. When the
Lagrangian decorrelation timescale is small, the turbulent flow takes less
time to forget its initial configuration compared to the typical time of
approach to the equilibrium of supersaturation. In the limit of
$\tau_{d}\rightarrow 0$, this argument suggests that $u(t)$ will become a
stochastic delta-correlated process, i.e., white noise. Indeed, by referring
to the theory of homogenisation [14, Chapter 11, Result 11.1], we are able to,
mathematically rigorously, reduce the two-dimensional system describing
$\mathrm{S}$ and $u$ to:

$\mathrm{d}\mathrm{S}=\left[-(B+C)\left(\mathrm{S}-\frac{C\mathrm{S}_{E}}{B+C}\right)\right](1+\mathrm{S})\mathrm{d}t+A\bar{u}(1+\mathrm{S})+A(1+\mathrm{S})\mathrm{d}W_{t},$
(17)

where a new constant $A$ has been introduced and defined as:

$A=a\sigma_{u}\sqrt{2\tau_{d}},$ (18)

which is the normalisation constant necessary to take this diffusion limit;
see [18, §6.3]. Because we have multiplicative noise, we have to specify the
stochastic calculus formalism being employed. For simplicity, Eq. (17) shall
be studied under the Itô formalism, although the Stratonovich viewpoint can be
taken instead [8]. For completeness, we note that in order to convert the
Stratonovich version of Eq. (17) into Itô, we apply the Itô-to-Stratonovich
correction $h(\mathrm{S})$ which reads as [8, 14]:

$h(\mathrm{S})=\frac{A^{2}}{2}\partial_{\mathrm{S}}(1+\mathrm{S})^{2}-\frac{A^{2}}{2}(1+\mathrm{S})\partial_{\mathrm{S}}(1+\mathrm{S})=\frac{A^{2}}{2}(1+\mathrm{S}).$
(19)

The resulting Itô stochastic differential equation is:

$\mathrm{d}\mathrm{S}=\left[-(B+C)\left(\mathrm{S}-\frac{2C\mathrm{S}_{E}-A^{2}}{2(B+C)}\right)\right](1+\mathrm{S})\mathrm{d}t+A\bar{u}(1+\mathrm{S})\mathrm{d}t+A(1+\mathrm{S})\mathrm{d}W_{t}.$
(20)

Notice that the equilibrium supersaturation $\mathrm{S}_{E}$ is now modified
by the term $-A^{2}$, resulting from the mere presence of multiplicative
noise. The best choice of stochastic formalism is not discussed here, although
this equation reveals that the full form of the Squire’s equation encodes
nonlinear interactions between supersaturation and turbulent fluctuations
which yield nontrivial corrections to the stochastic formulation of
supersaturation evolution.

Because the system is stochastic, the solutions of Eq. (17) give different
results for each noise realisation. To solve this problem, averages over all
possible realisations are taken so the system is described in terms of
probability distributions. In this section, we shall assume that in a
developing cloud the supersaturation budget quickly approaches stationarity,
so that the probability of encountering certain value of supersaturation is
provided the Fokker-Planck equation [17] associated with Eq. (17):

$\partial_{t}f(\mathrm{S},t)=\partial_{\mathrm{S}}\left[(B+C)\left(\mathrm{S}-\frac{C\mathrm{S}_{E}+A\bar{u}}{B+C}\right)(1+\mathrm{S})f(\mathrm{S},t)+\frac{A^{2}}{2}(1+\mathrm{S})^{2}\partial_{\mathrm{S}}f(\mathrm{S},t)\right],$
(21)

where $f(\mathrm{S},t)$ indicates— when normalised— the probability of
encountering a supersaturation of $\mathrm{S}$ at time $t$. Roughtly speaking,
this equation says that a density function is advected and diffused by the
linear and stochastic components of Eq. (17), respectively.

Unlike the analytical calculations of [7], the multiplicative noise and
nonlinearities involved in Eq. (17) suggest that the resulting probability
densities $f$ cannot be Gaussian. To find it out, we consider the stationary
version Eq. (21), which, after successive integrations detailed in Appendix B,
the normalised time-independent density is:

$\boxed{f_{1}(\mathrm{S})=\frac{\alpha^{\alpha+\beta+1}}{e^{\alpha}\Gamma(\alpha+\beta+1)}e^{-\alpha\mathrm{S}}(1+\mathrm{S})^{\alpha+\beta}},$
(22)

where we have introduced the nondimensional parameters $\alpha$, $\beta$ and
$\mathrm{S}^{\ast}$ defined as:

$\displaystyle\alpha$ $\displaystyle=\frac{2(B+C)}{A^{2}};$ (23a)
$\displaystyle\beta$ $\displaystyle=\alpha\mathrm{S}^{\ast};$ (23b)
$\displaystyle\mathrm{S}^{\ast}$
$\displaystyle=\frac{C\mathrm{S}_{E}+A\bar{u}}{B+C}.$ (23c)

The yielding moments can be computed accordingly in terms of successive Gamma
functions. The general formula for the uncentred moments of this distribution
is:

$\mathbb{E}\left[\mathrm{S}^{n}\right]=\sum_{k=0}^{n}(-1)^{k+1}\begin{pmatrix}n\\\
k\end{pmatrix}\alpha^{-k}\prod_{\ell=1}^{k}(\alpha+\beta+\ell).$ (24)

In particular, the formula for the variance is:

$\mathrm{Var}\left(\mathrm{S}\right)=\mathbb{E}\left[\left(\mathrm{S}-\mathbb{E}\left[\mathrm{S}\right]\right)^{2}\right]=\frac{\alpha+\beta+1}{\alpha^{2}}.$
(25)

The parameter $\alpha$ is non-dimensional and measures the relative strength
of the updraft fluctuations through $A$ and the combination of the absorption
and mixing timescales $B+C$. Hence, it follows that a large value of
$\alpha\gg 1$ yields a smaller variance. Moreover, in the limit of large
$\alpha$, the function $f_{1}$ can be recast into a Gaussian by means of
Laplace’s method; see, e.g., [4].

As $t$ tends to infinity, the expected value of $\mathrm{S}$ will converge to
an equilibrium value at a characteristic rate given by the phase relaxation
$\tau_{p}$ [13]. If the dynamics where deterministic, such rate would be given
by:

$\tau_{p}=\frac{1}{A\bar{u}+B+C},$ (26)

which is obtained by examining the exponent of the solution of the ordinary
differential equation (20); see details in [13]. In the stochastic context,
such relaxation rate is obtained by taking the expectation over all noise
realisations. However, it would be necessary to solve an open system involving
higher order moments of supersaturation [20, 23]. Hence, it is not clear
whether Eq. (26) is still the phase relaxation for noisy updraft fluctuations.
Nevertheless, the solution of the deterministic version of Eq. (20) is solved
and compared against an ensemble mean of Eq. (20), both cases subject to a
mean updraft of $\bar{u}=0.2\mathrm{m/s}$. Such ensemble mean is calculated by
taking $10^{4}$ noise realisations and averaging them at every time step. The
results are plotted in Figure 1, where we also include the relaxation curve
for the linear approximation that will be explained in the next section. The
initial condition of supersaturation is taken to be $0.5$ just for
demonstration purposes. It is observed that the deterministic relaxation
(orange curve) provides a good estimate for the stochastic one (blue curve),
although a slight divergence is seen as time increases.

####  3.1.1 Linear approximation

For stratocumulus and cumulus clouds, supersaturation does not go beyond
$0.5\%$, so that the approximation $1+\mathrm{S}\approx 1$ is widely taken in
the literature as a first order approximation of the Squire’s equation [10].
Under this approximation and some algebraic manipulations, Eq. (17) becomes an
OU process:

$\mathrm{d}\mathrm{S}=-(B+C)\left(\mathrm{S}-\mathrm{S}^{\ast}\right)\mathrm{d}t+A\mathrm{d}W_{t}.$
(27)

Such linear stochastic differential equation describes a process which has an
invariant density equal to a Gaussian distribution. In particular, mean and
variance are given by:

$\displaystyle\mathbb{E}\left[\mathrm{S}(t)\right]$
$\displaystyle=e^{-(B+C)t}\mathrm{S}(0)+\mathrm{S}^{\ast}\left(1-e^{-(B+C)t}\right)$
(28a) $\displaystyle\mathrm{Var}\left(\mathrm{S}(t)\right)$
$\displaystyle=\frac{1}{\alpha}\left(1-e^{-2(B+C)t}\right)=\frac{A^{2}}{2(B+C)}\left(1-e^{-2(B+C)t}\right)$
(28b)

As $t$ tends to infinity, the supersaturation value subject to a constant mean
updraft tend to an equilibrium value $\mathrm{S}^{\ast}$ exponentially fast,
with rate given by $-(B+C)$. Hence, the steady-state statistics are provided
by a Gaussian distribution:

$\boxed{f_{2}(\mathrm{S})=\frac{\alpha}{\sqrt{2\pi}}e^{-\frac{\alpha(\mathrm{S}-\mathrm{S}^{\ast})^{2}}{2}}}.$
(29)

Unlike the nonlinear supersaturation equation investigated in Section 3.1, the
phase relaxation for the present linear version yields the same relaxation
time as its deterministic analogue. This is given by:

$\tau_{p}=\frac{1}{B+C}.$ (30)

Note that the phase relaxation in Eq. (26) is different to Eq. (30), in that
the mean updraft $\bar{u}$ is not present in the latter. Consequently, it is
only in the limit of $\bar{u}\approx 0$ when both equations converge to
equilibrium at the same rate. For illustration, Eq. (28a) is plotted in Figure
1 in green colour, and demonstrates that the nonlinearity accelerates the
convergence to equilibrium.

Figure 1: Phase-relaxation to equilibrium. The blue curve represents the mean
of an ensemble of $10^{4}$ noise realisations of Eq. (20), with a prescribed
mean updraft of $\bar{u}=0.2$ $\mathrm{m/s}$. The orange curve is the solution
of Eq. (20), where noise is vanishing. The green curve is the exponential
relaxation time obtained from the linear equation (27). The ambient conditions
in this numerical experiment are: $\mathrm{S}(0)=0.5$, $T=-10^{\circ}$,
$p=50500\mathrm{Pa}$, $\bar{r}=5\mu\mathrm{m}$ and $N_{i}=100\mathrm{L}^{-1}$.

###  3.2 Slow decorrelation timescale: $\tau_{d}\gg 1$

When there is no timescale separation between the phase-relaxation and the
updraft fluctuations, the limit of $\tau_{d}$ going to zero cannot be taken
and therefore the white noise limit of the the previous section is not valid.
Because of this, the resulting nonlinear equation posses invariant statistics
that are intractable analytically. Hence, we start from Eq. (13) and assume,
as in Section 3.1.1, that $1+\mathrm{S}\approx 1$. We obtain the following
two-dimensional stochastic linear equation:

$\displaystyle\mathrm{d}\mathrm{S}$
$\displaystyle=-(B+C)\left(\mathrm{S}-\frac{C\mathrm{S}_{E}}{B+C}\right)\mathrm{d}t+au\mathrm{d}t,$
(31a) $\displaystyle\mathrm{d}u$
$\displaystyle=-\frac{1}{\tau_{d}}\left(u-\bar{u}\right)\mathrm{d}t+\sqrt{\frac{2}{\tau_{d}}}\sigma_{u}\mathrm{d}W_{t}.$
(31b)

Which can be compactly recast into matrix form:

$\mathrm{d}\begin{bmatrix}\mathrm{S}\\\
u\end{bmatrix}=\mathbf{B}\left(\begin{bmatrix}\mathrm{S}\\\
u\end{bmatrix}-\mathbf{m}\right)\mathrm{d}t+\mathbf{A}\mathrm{d}\mathbf{W}_{t},$
(32)

where $\mathbf{W}_{t}$ is a two-dimensional independent Wiener process and
where the matrices $\mathbf{B}$ and $\mathbf{A}$ and the vector $\mathbf{m}$
are defined as:

$\mathbf{B}=\begin{bmatrix}b_{11}&b_{12}\\\
0&b_{22}\end{bmatrix}=\begin{bmatrix}-(B+C)&a\\\ 0&-1/\tau_{d}\end{bmatrix},$
(33a) $\mathbf{A}=\begin{bmatrix}0&0\\\
0&a_{22}\end{bmatrix}=\begin{bmatrix}0&0\\\
0&(2/\tau_{d})^{1/2}\sigma_{u}\end{bmatrix},$ (33b)
$\mathbf{m}=\mathbf{B}^{-1}\begin{bmatrix}C\mathrm{S}_{E}\\\
\frac{\bar{u}}{\tau_{d}}\end{bmatrix}=\begin{bmatrix}\frac{C\mathrm{S}_{E}+a\bar{u}}{B+C}\\\
\bar{u}\end{bmatrix}.$ (33c)

Notice that the present equation is degenerate— noise affects directly to only
one variable—, although it will posses an invariant distribution wich
possesses a smooth density function. However, the two-dimensional covariance
matrix cannot be obtained straightforwardly out of the noise law, but by
computing the following matrix integral [14, Proposition 3.5]:

$\Sigma=\begin{bmatrix}\sigma_{11}&\sigma_{12}\\\
\sigma_{12}&\sigma_{22}\end{bmatrix}=\int_{0}^{\infty}e^{\mathbf{B}s}\mathbf{A}\mathbf{A}^{\top}e^{\mathbf{B}^{\top}s}\mathrm{d}s=\begin{bmatrix}\frac{a^{2}\sigma_{u}^{2}}{2(B+C)(B+C+1/\tau_{d})}&\frac{-a\sigma_{u}^{2}}{2(B+C+1/\tau_{d})}\\\
\frac{-a\sigma_{u}^{2}}{2(B+C+1/\tau_{d})}&\sigma_{u}^{2}\end{bmatrix}.$ (34)

The resulting process possesses, as in the one-dimensional case, a Gaussian
stationary distribution where now:

$\boxed{f_{3}(\mathrm{S},u)=(2\pi)^{-1}\left(\mathrm{det}\Sigma\right)^{-1/2}e^{-\frac{1}{2}\left([\mathrm{S},u]-\mathbf{m}^{\top}\right)\Sigma^{-1}\left([\mathrm{S},u]^{\top}-\mathbf{m}\right)}}.$
(35)

To derive the marginal distribution for $\mathrm{S}$ we apply the affine
transformation $P=[1,0]$ to the random variable $[\mathrm{S},u]^{\top}$ so
that $P[\mathrm{S},u]^{\top}=\mathrm{S}$. Hence,

$P\begin{bmatrix}\mathrm{S}\\\
u\end{bmatrix}=\mathrm{S}\sim\mathcal{N}\left(P\mathbf{m},P\Sigma
P^{\top}\right)=\mathcal{N}\left(\frac{C\mathrm{S}_{E}+a\bar{u}}{B+C},\frac{a^{2}\sigma_{u}^{2}}{2(B+C)(B+C+1/\tau_{d})}\right)$
(36)

The timescale separation argument to reduce the equation for the evolution of
supersaturation is only valid for small values of $\mathrm{TKE}$. This
hypothesis is tested and shown in Figure 2, where the variance of
supersaturation in a cold cloud is computed using the formulas here presented
and its numerical estimation using long time series of $10^{6}$ seconds. The
equations are integrated using a simple Euler-Maruyama method with a time-step
of $10^{-2}$ seconds [8].

Figure 2: Supersaturation variance as a function of TKE. The blue colours
correspond to the variances for the 2D linear system (31) with red-noise
updrafts, whereas the black colours correspond to the homogenised 1D equation
(27), where updrafts become white in time. The dots are obtained by
calculating the variance of $10^{6}$-second time series with a time-step of
$10^{-2}$ seconds, for each value of $\mathrm{TKE}$. The solid lines are the
predicted variance using the formulas (23a) and (36), for the blue and orange
curves, respectively. We highlight that when turbulence is less energetic,
both estimates become more similar. The ambient conditions in this numerical
experiment are: $\mathrm{S}_{E}=0$, $\bar{u}=0$, $T=-10^{\circ}$,
$p=50500\mathrm{Pa}$, $\bar{r}=5\mu\mathrm{m}$ and $N_{i}=100\mathrm{L}^{-1}$.

##  4 Diffusional growth and size fluctuations

In the Introduction we presented the full and closed equation for the
evolution of ice supersaturation accounting for the condensational growth of
droplets/ice crystals coupled to turbulent updraft fluctuations; see Eq. (7).
In Section 3, we investigated the properties of the stochastic Squires
equation, describing the evolution of supersaturation in the quasi-steady
case, where condensational growth is neglected. Here, we wish to explore the
effects of turbulent random forcings on the diffusional growth and the
latter’s on supersaturation evolution. Two default approximations are taken in
this section. First, we will assume that updraft fluctuations decorrelate
instantly so that the white-noise model is valid. Secondly, the small-
supersaturation approximation is taken: $1+\mathrm{S}\approx 1$. Then, by
applying the chain rule to Eq. (5), we can write the supersaturation equation
coupled to diffusional growth as:

$\displaystyle\mathrm{d}\mathrm{S}$
$\displaystyle=-C(\mathrm{S}-\mathrm{S}^{\ast}_{d})\mathrm{d}t-B_{d}r\mathrm{S}\mathrm{d}t+A\mathrm{d}W_{t}$
(37a) $\displaystyle\mathrm{d}r^{2}$
$\displaystyle=F\mathrm{S}\mathrm{d}t,\text{ }r^{2}>0$ (37b)
$\displaystyle\mathrm{d}r^{2}$ $\displaystyle=0,\text{ }r^{2}=0,\text{ and
}\mathrm{S}(t)<0.$ (37c)

Where some constants are introduced:

$\displaystyle\mathrm{S}^{\ast}_{d}$
$\displaystyle=\frac{C\mathrm{S}_{E}+A\bar{u}}{C};$ (38a) $\displaystyle
B_{d}$ $\displaystyle=bB_{0}N;$ (38b) $\displaystyle F$
$\displaystyle=2cA_{\circ}.$ (38c)

We note that the coupling of the equations (37a) and (37b) is nonlinear since
it involves the square-root of the variable $r^{2}$ and makes the analytics
intractable. However, we refer at this stage to Appendix C for the analysis of
the square-root stochastic process.

When the nonlinear term $r\mathrm{S}$ is small and close to zero, the mean and
variance of $r^{2}$ evolve according to:

$\displaystyle\mathbb{E}\left[r^{2}(t)\right]$
$\displaystyle=F\int^{t}\mathbb{E}\left[\mathrm{S}(s)\right]\mathrm{d}s+r^{2}(0)\approx\frac{Fe^{-Ct}}{C}\left(\mathrm{S}^{\ast}_{d}-\mathrm{S}(0)\right)+F\mathrm{S}^{\ast}_{d}t+r^{2}(0)$
(39a) $\displaystyle\mathbb{E}\left[\left(r^{2}(t)\right)^{2}\right]$
$\displaystyle=F^{2}\mathbb{E}\left[\int_{0}^{t}\int_{0}^{t}\mathrm{S}(s)\mathrm{S}(u)\mathrm{d}s\mathrm{d}u\right]=2F^{2}\int_{0}^{t}\int_{0}^{u}\frac{A^{2}}{2C}\left(e^{-C(u-s)}-e^{-C(u+s)}\right)\mathrm{d}s\mathrm{d}u$
(39b)
$\displaystyle=\frac{F^{2}A^{2}}{C^{2}}t-\frac{F^{2}A^{2}}{2C^{3}}+\frac{F^{2}A^{2}}{2C^{3}}\left(4e^{-Ct}-e^{-2Ct}\right).$
(39c)

These set of formulas are identical to those obtained in the 30s in the study
of Brownian motion— see [25]—, and imply that the evolution of the squared
radius follows Brownian paths where, in particular, the mean-squared
displacement scales linearly for large times, if negative values of squared
radii were allowed. Indeed, if $t\gg 1$ and the nonlinear coupling is small:

$\mathbb{E}\left[\left(r^{2}(t)\right)^{2}\right]\approx\frac{F^{2}A^{2}}{C^{2}}t.$
(40)

The boundary condition $r^{2}=0$ is strictly necessary since it is possible
that trajectories of $r(t)$ in Eq. (37b) vanish. When that happens, it means
that the particles or droplets in question have evaporated and that the
formula for the variance in Eq. (39c) has to be reinitialised once the
droplets and particles have reactivated.

While under this framework there is no stationary distribution with finite
variance for particle radius, it was shown in [23] that if the ambient
supersaturation $\mathrm{S}_{E}$ is negative, i.e., subsaturated, the
probability distribution of $r^{2}$ will possess the structure of an
exponential function with an Dirac-peak located at $r^{2}=0$ which arises from
the boundary condition in Eq. (37c).

###  4.1 Fluctuations in droplet radius

In the previous section we clarified that if supersaturation is let to be
driven by random turbulent updrafts, the mean-square radius grows linearly in
time and, therefore, unbounded Brownian excursions can be expected when
solving the condensational growth equation. When the cloud in question is in
contact with a subsaturated environment or the system allows for droplet
evaporation or sedimentation, it is expected that the trajectories in the
$\mathrm{S}$-$r^{2}$-plane of Eq. (37) will display cycles, where $r^{2}$
grows but then vanishes and sticks at the boundary of $r^{2}=0$ for an open
interval of time; see [23, Figure 4]. In this section, we aim at calculating
the net effects of the mentioned cycles, by modelling them as an extra random
forcing at the microphysical term. For this purpose we introduce a new
constant $\sigma_{r}$ which indicates the standard deviation of droplet radius
fluctuations. For small particles where condensation is a dominating growth
factor, $\sigma_{r}$ must be proportional to $F$ and the standard deviation of
$\mathrm{S}$. To simplify the expression, we shall present the results for
$\mathrm{S}_{E}=0$ and $\bar{u}=0$.

The starting point is the $1+\mathrm{S}\approx 1$ approximation that together
with fluctuations in droplet size lead to a new model for supersaturation:

$\mathrm{d}\mathrm{S}(t)=-C\mathrm{S}(t)\mathrm{d}t-B_{d}\mathrm{S}(t)\left(\bar{r}\mathrm{d}t+\sigma_{r}\mathrm{d}W^{(1)}_{t}\right)+A\mathrm{d}W^{(2)}_{t},$
(41)

where $W^{(1)}$ and $W^{(2)}$ are two independent Wiener processes. Also, the
parameter $B$ has been modified to $B=bB_{0}N$. In this case, solving the
stationary Fokker-Planck equation for $\mathrm{S}^{\ast}=0$ yields the
following non-normalised stationary distribution:

$\boxed{f_{4}(\mathrm{S})=\left(A^{2}+\sigma_{r}^{2}B_{d}^{2}\mathrm{S}^{2}\right)^{-1-\frac{C+B_{d}\bar{r}}{B_{d}^{2}\sigma_{r}^{2}}}}.$
(42)

This formula constitutes a power-law, for which higher momenta might not be
well defined. However, by examining the exponent, it follows that the
distribution $f_{4}$ has $M$ momenta, if
$M/2<1+(C+B_{d}\bar{r})/B_{d}^{2}\sigma_{r}^{2}$. The function $f_{4}$
describes the steady-state statistics when $W^{(1)}_{t}\neq W^{(2)}_{t}$, but
since fluctuations in droplets size typically originate in fluctuations of
supersaturation, it is important to consider the the case with
$W^{(1)}_{t}=W^{(2)}_{t}$, so that noise sources are correlated. In such
scenario, the stationary distribution is:

$\boxed{f_{5}(\mathrm{S})=e^{-\frac{2(C+B_{d}\bar{r})A}{B_{d}^{2}\sigma_{r}^{2}}(A-B_{d}\sigma_{r}\mathrm{S})^{-1}}\left(A-B_{d}\sigma_{r}\mathrm{S}\right)^{-2-\frac{2(C+B_{d}\bar{r})}{B_{d}^{2}\sigma_{r}^{2}}}}.$
(43)

The probability distributions of Eq. (42) and (43) differ from a Gaussian
distribution since the noise appears in a multiplicative way. It is clear,
though, that the stochastic process Eq. (41) will converge to linear evolution
of supersaturation as $\sigma_{r}$ tends to zero, which yields the Gaussian
distribution of Eq. (29). It is possible, on the other hand, to show the
pointwise convergence of Eq. (42) and (43) to Eq. (29), by means of Laplace’s
method; see, e.g., [4].

##  5 Parametrisation of mixed-phase clouds

The relevance of the analytical calculations done in the previous sections is
discussed here in the context of mixed-phase clouds. The coexistence of liquid
water and ice is thermodynamically unstable so that, in freezing temperatures,
ice will inevitably grow at the expense of liquid water [3]. Observations, on
the other hand, indicate that mixed-phase conditions are not just a transient
microphysical state, but that supercooled liquid water can be maintained at
cloud top temperatures down to $-40^{\circ}$. The activation of supercooled
liquid water in a cloud parcel is in general achieved if the following two
criteria are met: (i) the vertical velocity must exceed a threshold value and
(ii) the cloud parcel must be lifted to a threshold altitude. Under this
context, mixed-phased conditions can be kept for a long term [12]. When
considering an icy cloud, Eq. (13) has to be interpreted with respect to ice
supersaturation $\mathrm{S}_{i}$.

Under this theoretical setting, [7] proposed a method to estimate the variance
of ice supersaturation using a Gaussian probability distribution, here
revisited in Section 3.1.1. Indeed, when supersaturation with respect to ice
exceeds the value at liquid water saturation, liquid water is activated.
Supersaturation with respect to ice relates to that of water,
$\mathrm{S}_{w}$, in terms of the ratio of their respective equilibrium vapour
pressures:

$\mathrm{S}_{i}=\eta(T)\mathrm{S}_{w}+\eta(T)-1.$ (44)

where $\eta(T)=E_{w}(T)/E_{i}(T)$. As a consequence, supersaturation with
respect to ice at water saturation $\mathrm{S}_{iw}$ is given by:

$\mathrm{S}_{iw}=\eta(T)-1.$ (45)

Hence, the fraction of cloud in mixed-phase conditions and liquid water
content are determined by the tails of the probability density functions in
$\\{f_{k}\\}_{k=1}^{5}$ from $\mathrm{S}_{iw}$ to infinity. Concretely, the
fraction of cloud parcel that is taken to be in mixed-phase is the total time
$\mathrm{S}_{i}$ spends above $\mathrm{S}_{iw}$ or, in other words— by
invoking ergodicity—, the integral of $f_{k}$ from $\mathrm{S}_{iw}$ to
infinity:

$\mathrm{C}_{f_{k}}=\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{t}\mathbf{1}_{\mathrm{S}_{i}\geq\mathrm{S}_{iw}}\left(\mathrm{S}_{i}(s)\right)\mathrm{d}s=\int_{\mathrm{S}_{iw}}^{\infty}f_{k}(\mathrm{S}_{i})\mathrm{d}\mathrm{S}_{i},$
(46)

where $\mathbf{1}_{\mathrm{S}_{i}\geq\mathrm{S}_{iw}}$ is the characteristic
function for values of $\mathrm{S}_{i}$ larger than $\mathrm{S}_{iw}$. As
noted in [7], in the Gaussian case the cloud fraction is given explicitly by:

$\mathrm{C}_{f_{2}}=\frac{1}{2}\mathrm{erfc}\left(\sqrt{\frac{\alpha}{2}}\left(\mathrm{S}_{iw}-\mathrm{S}^{\ast}\right)\right).$
(47)

Moreover, the parcel’s LWC is calculated by assuming that ice supersaturation
above that with respect to liquid water is converted into droplets. The LWC
is, then, estimated by integrating the adjustment formula of Eq. (44) to
obtain the condensed water amount:

$\left\langle
q\right\rangle_{k}=\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{t}\mathbf{1}_{\mathrm{S}_{i}\geq\mathrm{S}_{iw}}(\mathrm{S}_{i}(s))\left(\mathrm{S}_{i}(s)-\mathrm{S}_{iw}\right)\rho_{a}q_{si}\mathrm{d}\mathrm{S}_{i}=\int_{\mathrm{S}_{iw}}^{\infty}\left(\mathrm{S}_{i}-\mathrm{S}_{iw}\right)f_{k}(\mathrm{S}_{i})\rho_{a}q_{si}\mathrm{d}\mathrm{S}_{i}$
(48)

where $\rho_{a}q_{si}$ is the factor that converts liquid water content to
supersaturation values.

The five probability distributions obtained in the previous sections are now
used to compute the cloud fraction and the mean LWC of a subgrid cloud parcel
for a range of free parameters. Such free parameters are the $\mathrm{TKE}$,
as a proxy for turbulent forcing, and variance of the droplet radius
fluctuations. In Figure 3(a) and (b) we show the dependence of the mentioned
partial moments on values of $\mathrm{TKE}$ at two temperatures indicated in
the captions. Such statistics where computed using a simple quadrature scheme
on the interval $[-10,10]$ so that all the considered PDFs integrate to unity
with a tolerance of $10^{-10}$. We observe a monotone dependence on
$\mathrm{TKE}$ in all the PDFs but for the Gamma distribution, which yielded
decreasing cloud fractions for values of $\mathrm{TKE}>6\mathrm{m^{2}s^{-2}}$
at $-10^{\circ}$ and $\mathrm{TKE}>4\mathrm{m^{2}s^{-2}}$ at $-5^{\circ}$.
Such change in trend is due to the displacement of the mode and tail-thickness
in the Gamma distribution as the location factor is altered due to the
multiplicative noise. The nonlinear interaction of updraft fluctuations and
supersaturation, hence, prevent the increase in mixed-phase conditions for
large values of turbulent forcing strength.

When droplet size fluctuations are allowed, the first three calculated
distributions $f_{1},f_{2}$ and $f_{3}$ naturally yield the same statistics
for cloud fraction and mean LWC. Contrarily, when noisy variations in radius
are allowed the distributions of $f_{4}$ and $f_{5}$ are likely to display a
dependence on $\sigma_{r}$. This dependence is shown in Figure 4, where the
cloud fraction and LWC are calculated as a function $\sigma_{r}$, for
$f_{2},f_{4}$ and $f_{5}$. Because $f_{4}$ and $f_{5}$ are expensive to
evaluate at small values of $\sigma_{r}$, ergodic averages are computed
instead, i.e., the first equalities in Eq. (46) and (48). Such averages are
taken over an integration of Eq. (41) over $10^{6}$ seconds. While $f_{2}$ is
expectedly constant, $f_{4}$ and $f_{5}$ start to be dependent for values
larger than $\sigma^{2}_{r}\geq 10^{-12}\mathrm{m}^{2}$, which corresponds to
a standard deviation $20\%$ of the mean droplet radius $\bar{r}$. In
particular, it is observed that $\mathrm{C}_{f_{4}}$ grows up to $1\%$ when
droplet fluctuations have a variance of $10^{-8}\mathrm{m}$. On the other
hand, $\mathrm{C}_{f_{5}}$ appears to be independent of $\sigma_{r}$. We
recall that this case corresponds to when the fluctuations in $r$ are
independent of the noisy updraft. Such independence does not hold in LWC,
where, $\langle q\rangle_{5}$ appears to grow, as a consequence of the
fattening of the tails.

(a)

(b)

Figure 3: Cloud fraction and mean LWC as a function of TKE. The cloud fraction
and mean LWC are calculated using Eq. (46) and Eq. (48), respectively, against
$\mathrm{TKE}$, for each analytical supersaturation distribution
$\\{f_{k}\\}_{k=1}^{5}$ and for two temperature configurations. In (a), a
temperature of $-10^{\circ}$ is considered, in (b), $-5^{\circ}$. The rest of
the ambient conditions are: $\mathrm{S}_{E}=0$, $\bar{u}=0$,
$\sigma^{2}_{r}=10^{-6}\mathrm{m}^{2}$, $p=50500\mathrm{Pa}$,
$\bar{r}=5\mu\mathrm{m}$ and $N_{i}=100\mathrm{L}^{-1}$.

(a)

(b)

Figure 4: Cloud fraction and mean LWC as a function of $\sigma_{r}$. The cloud
fraction and mean LWC are calculated using Eq. (46) and Eq. (48),
respectively, against $\sigma_{r}$, for the analytical distributions
$f_{1},f_{4}$ and $f_{5}$ and for two temperature configurations. For small
values of $\sigma_{r}$, the functions $f_{4}$ and $f_{5}$ are computationally
expensive to evaluate so ergodic averages of $10^{6}$ seconds are taken
instead. In (a), a temperature of $-10^{\circ}$ is considered, in (b),
$-5^{\circ}$. The rest of the ambient conditions are: $\mathrm{S}_{E}=0$,
$\bar{u}=0$, $\mathrm{TKE}=2.612\mathrm{m^{2}/s^{2}}$, $p=50500\mathrm{Pa}$,
$\bar{r}=5\mu\mathrm{m}$ and $N_{i}=100\mathrm{L}^{-1}$.

##  6 Discussion

In this paper, the analysis of the stochastic Squires equation has been done,
with the aim of providing analytical formulas for the distribution of
supersaturation in a cloud parcel. Such equation describes the evolution of
supersaturation over time by, essentially, taking into account the sources and
sinks of relative humidity due to (a), adiabatic cooling via turbulent
updrafts and (b), water vapour condensation onto droplets/particles. The
former are here modelled as a red-noise process, in accordance with the
spatial and time correlations of isotropic fluctuations in the inertial regime
[18]. Secondly, droplet growth by vapour condensation is studied and noted— as
firstly done in [23]— that the variance of squared droplet radius grows
linearly in time and, hence, it does not have bounded steady-state statistics.
Here, we propose to include condensational growth effects by adding an
explicit stochastic term for droplet size fluctuations which yield
analytically tractable probability distributions for supersaturation. In
total, a number of five different probability distributions where computed.

In the present paper, the theoretical results of [7]— here condensed in
Section 3.1.1— have been generalised to a wider range of contexts. First of
all, the Squires equation is considered in Section 2 in its full nonlinear
version and has been shown to possess a Gamma-like stationary distribution—
here denoted as $f_{1}$— that deviates from a Gaussian according to the
parameter $\alpha$. Such parameter, also found in [7], is a non-dimensional
ratio between the strength of turbulent fluctuations and the phase relaxation
coefficient. Thus, as turbulent fluctuations decrease in variance (relative to
the microphysical or mixing timescale), $f_{1}$ becomes better and better
approximated by the Gaussian distribution $f_{2}$. Indeed, that can be seen
from the computation of partial moments— cloud fraction and mean LWC— in
Figure 3.

The timescale separation assumption, needed to compute $f_{1}$ and $f_{2}$, is
lifted in Section 3.2. Indeed, $f_{2}$ is only valid when times are much
greater that the Lagrangian decorrelation timescale so that updraft
fluctuations become white in time. The distribution $f_{3}$ is still Gaussian,
but it does not invoke such timescale separation between updraft fluctuations
and supersaturation. On Figure 2, we show that the variance predicted by
$f_{2}$ and $f_{3}$ diverge as $\mathrm{TKE}$ increases.

The main assumption needed to compute $f_{1},f_{2}$ and $f_{3}$ is the quasi-
steady approximation, whereby the droplet or ice particle radius is considered
constant. It was shown in [23] that the long-term variance of droplet squared-
radius scales linearly in time, similar to Brownian motion. This is revisited
here in Section 4. We argue that such result is only valid for short times,
since large families of droplets are subject to processes like sedimentation,
evaporation or mixing with exterior dry air that provoke a _memory loss_ in
collective droplet growth and, hence, the variance is reinitialised. In this
work, instead, we investigate random fluctuations in droplet size in two
contexts: (i) fluctuations in updrafts are uncorrelated to those of droplet
size and (ii), the source of noise is the same, albeit with different
intensities. The net effect of droplet radius fluctuations in summarised in
Figure 4, where the cloud fraction and mean LWC are computed as a function of
$\sigma_{r}$. It is found that correlated radius and updraft fluctuations
yield a more sensitive probability distribution, which deviates severely from
the quasi-steady approximation by up to $+5\%$ in cloud fraction when
$\sigma_{r}\geq 10^{-6}\mathrm{m}$. Contrarily, the mean LWC is negatively
correlated with fluctuations in $\sigma_{r}$ in case of $f_{4}$. Regarding the
uncorrelated sources of noise, the cloud fraction appears to be weakly
dependant in $\sigma_{r}$. On the other hand, mean LWC correlates positively
is droplet size fluctuation variance.

The present stochastic analysis of supersaturation is argued to be a useful
framework to find probabilistic formulas for the parametrisation of subgrid
cloud properties. However, a deeper investigation of Eq. (7) would be useful
in this stochastic framework. One step forward would be to impose a
characteristic time for the loss of memory, so that the integro-differential
equation can be replaced by a simpler expression, possibly some form of noise
with suitable time-decorrelation properties. We anticipate that this would
entail technical difficulties due to the square-root nonlinearity— here
minimally tackled in Appendix C—, so research should be oriented in this
direction. In general, we belief that this approach can be extended to more
general contexts, possibly, by including more microphysical processes that
affect the growth of liquid droplets or ice particles and, hence, the overall
regulation of a cloud supersaturation budget.

## Data availability

The data to produce the figures in this paper is obtained from the numerical
integration of the equations in question. The corresponding code is available
upon request.

## Acknowledgements

The authors would like to thank S. Roncoroni, P. Field and B. Devenish for
their comments, suggestions and kind reception at the MetOffice UK. MSG is
grateful to the Mathematics of Planet Earth Centre for Doctoral Training (MPE
CDT) for making this collaboration possible. MSG acknowledges and is grateful
for the support of the Institute of Mathematics and its Applications (grant
number: SGS21/08). MSG is thankful to I. Koren, M. D. Chekroun, the cloud
physics group and the graduate school at the Weizmann Institute of Science for
providing a most inspiring environment.

##  Appendix A List of some used notations and symbols

Symbol | Units | Description
---|---|---
$A$ | [$\mathrm{s}^{-1/2}$] | $a_{i}\sigma_{u}\tau_{d}^{1/2}$
$A_{i}$ | [$\mathrm{m}^{2}\mathrm{s^{-1}}$] | $\left(\frac{\rho_{i}L_{i}^{2}}{kR^{v}T^{2}}+\frac{\rho_{i}R_{v}T}{E_{i}D}\right)^{-1}$
$A_{w}$ | [$\mathrm{m}^{2}\mathrm{s^{-1}}$] | $\left(\frac{\rho_{w}L_{w}^{2}}{kR^{v}T^{2}}+\frac{\rho_{w}R_{v}T}{E_{w}D}\right)^{-1}$
$A_{\circ}$ | [$\mathrm{m}^{2}\mathrm{s^{-1}}$] |  $A_{i}$ or $A_{w}$ depending on the context
$a_{w}$ | [$\mathrm{m}^{-1}$] | $\frac{g}{R_{a}T}\left(\frac{L_{w}R_{a}}{c_{p}R_{v}T}-1\right)$
$a_{i}$ | [$\mathrm{m}^{-1}$] | $\frac{g}{R_{a}T}\left(\frac{L_{i}R_{a}}{c_{p}R_{v}T}-1\right)$
$a$ | [$\mathrm{m}^{-1}$] |  $a_{i}$ or $a_{w}$ depending on the context
$B$ | $[\mathrm{s}^{-1}]$ | $b_{i}B_{0}N\bar{r}$
$B_{d}$ | $[\mathrm{r^{-1}~{}s}^{-1}]$ | $b_{i}B_{0}N$
$b_{w}$ | $[-]$ | $\frac{1}{q_{v}}+\frac{L_{w}^{2}}{c_{p}R_{v}T^{2}}$
$b_{i}$ | $[-]$ | $\frac{1}{q_{v}}+\frac{L_{i}^{2}}{c_{p}R_{v}T^{2}}$
$b_{m}$ | $[-]$ | $\frac{1}{q_{v}}+\frac{L_{i}L_{w}}{c_{p}R_{v}T^{2}}$
$b$ | $[-]$ |  $b_{w}$, $b_{i}$ or $b_{m}$ depending on the context
$B_{0}$ | $[\mathrm{m}^{2}\mathrm{s^{-1}}]$ | $\frac{4\pi\rho_{i}A_{i}c}{\rho_{a}}$
$C_{0}$ | $[-]$ | Lagrangian structure function constant($=10$)
$C_{f}$ | $[-]$ | Mixed-phase cloud fraction
$c$ | $[-]$ | Ice-crystal capacitance, shape factor ($=1$)
$c_{p}$ | $[\mathrm{J}\mathrm{kg}^{-1}\mathrm{K}^{-1}]$ | Specific heat capacity of moist air at constant pressure
$D$ | $[\mathrm{m}^{2}\mathrm{s^{-1}}]$ | Water-vapour diffusion coefficient in air
$e$ | $[\mathrm{Pa}]$ | water vapour pressure
Table 1: List of symbols. Symbol | Units | Description
---|---|---
$E_{i}$ | $[\mathrm{Pa}]$ | Saturation vapour pressure over ice
$E_{w}$ | $[\mathrm{Pa}]$ | Saturation vapour pressure over liquid water
$E$ | $[\mathrm{Pa}]$ | Saturation vapour pressure over liquid or ice depending on the context
$\varepsilon$ | $[\mathrm{m}^{2}\mathrm{s}^{-3}]$ | Eddy dissipation rate
$F$ | $[\mathrm{m}^{2}\mathrm{s}^{-1}]$ | $2cA_{\circ}$
$g$ | $[\mathrm{m}\mathrm{s^{-2}}]$ | Acceleration due to gravity
$k$ | $[\mathrm{J}\mathrm{m}^{-1}\mathrm{s}^{-1}\mathrm{K}^{-1}]$ | Heat conductivity coefficient in air
$L$ | $[\mathrm{m}]$ | Vertical length of turbulent zone
$L_{i}$ | $[\mathrm{J}\mathrm{kg}^{-1}]$ | Latent heat for sublimation plus melting
$L_{w}$ | $[\mathrm{J}\mathrm{kg}^{-1}]$ | Latent heat for vaporisation
$N_{i}$ | $[\mathrm{L}^{-1}]$ | Concentration of ice particles
$N_{w}$ | $[\mathrm{L}^{-1}]$ | Concentration of water droplets
$N$ | $[\mathrm{L}^{-1}]$ | Concentration of water droplets/ice particles depending on the context
$p$ | $[\mathrm{Pa}]$ | Pressure
$\langle q\rangle$ | $[\mathrm{kg~{}m}^{-3}]$ | Domain mean liquid water content
$q_{si}$ | $[\mathrm{kg~{}kg}^{-1}]$ | Mass mixing ratio of vapour at ice saturation
$q_{v}$ | $[\mathrm{kg~{}kg^{-1}}]$ | Mass mixing ratio of water vapour
$q$ | $[\mathrm{kg~{}kg^{-1}}]$ | Liquid water or ice mixing ratio depending on the context
$R_{v},R_{a}$ | $[\mathrm{J~{}kg^{-1}~{}K^{-1}}]$ | Specific gas constant of water vapour and air, respectively
Table 2: List of symbols. Symbol | Units | Description
---|---|---
$\mathrm{S}_{iw}$ | $[-]$ | $E_{w}/E_{i}-1$
$\mathrm{S}_{w}$ | $[-]$ | Supersaturation with respect to liquid water
$\mathrm{S}_{i}$ | $[-]$ | Supersaturation with respect to ice
$\mathrm{S}$ | $[-]$ | Supersaturation with respect to liquid water or ice depending on the context
$\mathrm{S}_{E}$ | $[-]$ | Environmental supersaturation with respect to liquid water or ice
$\mathrm{S}_{i}^{\ast}$ | $[-]$ | $\mathrm{S}_{iw}+q/(\rho_{a}q_{si})$
$T$ | $[\mathrm{K}]$ | Temperature
$\mathrm{TKE}$ | $[\mathrm{m^{2}/s^{2}}]$ | Turbulent kinetic energy
$u$ | $[\mathrm{m~{}s^{-1}}]$ | Vertical velocity
$\bar{u}$ | $[\mathrm{m~{}s^{-1}}]$ | Mean vertical velocity
$\epsilon$ | $[\mathrm{m^{2}s^{-3}}]$ | Eddy dissipation rate
$\nu$ | $[\mathrm{m^{2}~{}s^{-1}}]$ | Kinematic viscosity
$\rho_{i}$ | $[\mathrm{kg~{}m^{-3}}]$ | Effective density of ice
$\rho_{a}$ | $[\mathrm{kg~{}m^{-3}}]$ | Density of air
$r$ | $[\mathrm{m}]$ | Droplet or ice particle radius depending on the context
$\bar{r}$ | $[\mathrm{m}]$ | Mean droplet or ice particle radius depending on the context
$\sigma_{u}$ | $[\mathrm{m~{}s^{-1}}]$ | Standard deviation of vertical velocity fluctuations
$\sigma_{r}$ | $[\mathrm{m}]$ | Standard deviation of droplet/particle size fluctuations
$\tau_{d}$ | $[\mathrm{s}]$ | Lagrangian decorrelation timescale
Table 3: List of symbols.

##  Appendix B Stationary supersaturation distribution

The calculation of the stationary distributions of each case study is done by
studying the Fokker-Planck representation of the stochastic processes in
question [17]. Such equation describes how probability distributions evolve in
time towards its stationary state. For simplicity, in this appendix we will
just show how to derive, step by step, the distribution for ice
supersaturation Eq. (22). In this case, the Fokker-Planck equation associated
with Eq. (17) is given by Eq. (21). Then, its stationary version, for
$\mathrm{S}_{E}=\bar{u}=0$, reads as:

$0=\partial_{\mathrm{S}}\left[(B+C)\mathrm{S}(1+\mathrm{S})f(\mathrm{S})+\partial_{\mathrm{S}}\left(\frac{A^{2}}{2}(1+\mathrm{S})^{2}f(\mathrm{S})\right)\right],$
(B.1)

where the time-dependence has been dropped. Integrating from $\ell$ (yet to be
determined) to $\mathrm{S}$ and assuming that the stationary distribution $f$
and $\partial_{\mathrm{S}}f$ vanish at $\ell$:

$\displaystyle
0=(B+C)\mathrm{S}(1+\mathrm{S})f(\mathrm{S})+A^{2}(1+\mathrm{S})f(\mathrm{S})+\frac{A^{2}}{2}(1+\mathrm{S})^{2}\partial_{\mathrm{S}}f(\mathrm{S}).$
(B.2)

We now rearrange the equation to make it homogeneous on both sides:

$\frac{\partial_{\mathrm{S}}f(\mathrm{S})}{f(\mathrm{S})}=-\frac{2(B+C)\mathrm{S}}{A^{2}(1+\mathrm{S})}-\frac{2}{1+\mathrm{S}}.$
(B.3)

Integrating on both sides:

$\mathrm{log}(f(\mathrm{S}))=\frac{2(B+C)}{A^{2}}\mathrm{log}(1+\mathrm{S})-\frac{2(B+C)}{A^{2}}\mathrm{S}-2\mathrm{log}(1+\mathrm{S}).$
(B.4)

Taking exponentials on both sides:

$f(\mathrm{S})=e^{-\frac{2(B+C)}{A^{2}}\mathrm{S}}(1+\mathrm{S})^{-2+\frac{2(B+C)}{A^{2}}}.$
(B.5)

This is a non-normalised solution for the stationary Fokker-Planck equation,
which has a singularity at $-1$. As $\mathrm{S}$ tends to infinity,
$f(\mathrm{S})$ goes to zero asymptotically. We now calculate the
normalisation constant $C$, for which we employ the parameter $\alpha$ like in
Eq. (23a):

$\displaystyle
C=\int_{-1}^{\infty}e^{\alpha(1+\mathrm{S})}(1+\mathrm{S})^{-2+\alpha}\mathrm{d}\mathrm{S}=e^{\alpha}\alpha^{1-\alpha}\int_{0}^{\infty}e^{-z}z^{\alpha-2}\mathrm{d}z=e^{-\alpha}\alpha^{1-\alpha}\Gamma(\alpha-1).$
(B.6)

Finally, $f/C$ gives Eq. (22). This process is repeated for every probability
distribution $\\{f_{k}\\}_{k=1}^{5}$, although the general formula is derived
below.

###  B.1 The general case

A general one-dimensional SDE reads as:

$\mathrm{d}x=V(x)\mathrm{d}t+\sigma(x)\mathrm{d}W_{t},$ (B.7)

where $W_{t}$ is a standard Wiener process and $V$ and $\sigma\neq 0$ have the
regularity so that solutions distribute according to a smooth probability
density function [8, 14]. The associated Fokker-Planck equation is:

$\partial_{t}f=\partial_{x}\left[-V(x)f+\frac{1}{2}\partial_{x}\left(\sigma^{2}(x)f\right)\right],$
(B.8)

where $f$ is a probability density and a function of $x$ and $t$. The
stationary distribution of Eq. (B.7) is obtained by setting $\partial_{t}f=0$,
and solving for $f$:

$\displaystyle
0=\partial_{x}\left[-V(x)f+\frac{1}{2}\partial_{x}\left(\sigma^{2}(x)f\right)\right].$
(B.9)

Integrating from $\ell$, where $f$ is assumed to vanish, to $x$ and
rearranging to make the equation homogeneous,

$\frac{\partial_{x}f(x)}{f(x)}=\frac{2V(x)}{\sigma^{2}(x)}-\frac{\partial_{x}\left(\sigma^{2}(x)\right)}{\sigma^{2}(x)}$
(B.10)

solving the indefinite integral on both sides yields:

$\log
f(x)=\int^{x}\frac{2V(x)}{\sigma^{2}(x)}\mathrm{d}x-\log\left(\sigma^{2}(x)\right),$
(B.11)

where “$\int^{x}$” denotes the indefinite integral. We now take the
exponentials on both sides:

$f(x)=e^{\int^{x}\frac{2V(x)}{\sigma^{2}(x)}}\mathrm{d}x\sigma^{-2}(x).$
(B.12)

One is left with finding the normalising constant.

##  Appendix C The square-root process

The condensational growth equation is easily integrated into the
supersaturation equation, although its expression depends on the square-root
of the initial particle size plus its fluctuating part due to supersaturation
variations; see Eq. (7). If such fluctuations, here denoted as
$\delta_{\mathrm{S}}(t)$, decorrelate instantly and have zero mean, such
expression is rewritten as:

$r(t)=\sqrt{r^{2}(t)}=\sqrt{r(0)^{2}+2cA_{\circ}\int_{0}^{t}\mathrm{S}(s)\mathrm{d}s}\approx\sqrt{r(0)^{2}+\delta_{\mathrm{S}}(t)},$
(C.1)

where the variance of $\delta_{\mathrm{S}}(t)$, $\sigma^{2}_{r^{2}}$, is
proportional to the diffusional constant $A_{\circ}$. Deriving the statistics
of a squared-root stochastic process is difficult, although in the limit of
variance fluctuations of $r^{2}$ begin small, or when the diffusional constant
$A_{\circ}$ is small, we can derive the steady-state mean and variance of the
square-root process as an expansion:

$\displaystyle\mathbb{E}\left[r\right]$
$\displaystyle=\mathbb{E}\left[\sqrt{r^{2}}\right]=r(0)-\frac{1}{8}r(0)^{-3}\sigma^{2}_{r^{2}}+\mathcal{O}\left(\sigma_{r^{2}}^{4}\right)$
(C.2a) $\displaystyle\mathrm{Var}\left[r\right]$
$\displaystyle=\mathrm{Var}\left[\sqrt{r^{2}}\right]=\mathbb{E}\left[r^{2}\right]-\mathbb{E}\left[\sqrt{r^{2}}\right]^{2}=\frac{1}{4r(0)^{2}}\sigma^{2}_{r^{2}}+\mathcal{O}\left(\sigma^{4}_{r^{2}}\right).$
(C.2b)

Surprisingly, the variance of the square-root process only scales inversely
proportionally to $r(0)^{2}$. This approach assumes that $r(0)$ is the mean
radius and that fluctuations $\delta_{\mathrm{S}}(t)$ are so small that
$r^{2}$ remains positive. This result is general and can be applied to any
square-root random variable, with positive mean, and in the limit of small
variance.

To support this analytical expansion, we numerically sampled an adimensional
random variable $X$, normally distributed with mean $2$ and standard-deviation
$\sigma_{X}$, where the latter takes 250 equispaced values between $10^{-4}$
and $0.5$. The sample is of size $10^{5}$ draws. With this set of data, we are
able to numerically estimate the mean and variance of the square root random
variable, $\sqrt{X}$, discarding all samples that gave negative values. The
aim is to predict the variance and mean of $\sqrt{X}$, using the truncated
expansions of Eq. (C.2) and the prescribed values of $\sigma_{X}$. The
truncation is done at $\mathcal{O}\left(\sigma^{4}_{X}\right)$. Indeed, in
Figure 5 the analytical predictions— plotted in solid coloured curves— match
to a high degree of accuracy the numerically sampled random variable
$\sqrt{X}$— plotted in black dots—. When $\sigma_{X}$ becomes larger, a
moderate deviation is observed, as expected.

Figure 5: Variance and mean of the square-root random variable. The black dots
are calculated as follows: for each value of $\sigma_{X}$, the variance and
mean of the random variable $\sqrt{X}$ are numerically estimated by taking the
square root of $10^{5}$ draws of a normal random variable $X$ with mean $2$
and standard-deviation $\sigma_{X}$. The blue and orange solid lines indicate
the truncated predictions of Eq. (C.2) for the variance and mean,
respectively.

## References

  * [1] A. Arakawa and W. H. Schubert.  Interaction of a cumulus cloud ensemble with the large-scale environment, part i.  Journal of Atmospheric Sciences, 31(3):674 – 701, 1974.

  * [2] J. T. Bartlett and P. R. Jonas.  On the dispersion of the sizes of droplets growing by condensation in turbulent clouds.  Quarterly Journal of the Royal Meteorological Society, 98(415):150–164, 1972.
  * [3] T. Bergeron.  On the physics of clouds and precipitation, pages 156––178.  Proces Verbaux de l’Association de Meteorologie. International Union of Geodesy and Geophysics, 1935.
  * [4] R. W. Butler.  Exponential families and tilted distributions, page 145–182.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2007.
  * [5] B. J. Devenish, K. Furtado, and D. J. Thomson.  Analytical solutions of the supersaturation equation for a warm cloud.  Journal of the Atmospheric Sciences, 73(9):3453–3465, 2016.
  * [6] E. Eytan, A. Khain, M. Pinsky, O. Altaratz, J. Shpund, and I. Koren.  Shallow cumulus properties as captured by adiabatic fraction in high-resolution les simulations.  Journal of the Atmospheric Sciences, 79(2):409 – 428, 2022.
  * [7] P. R. Field, A. A. Hill, K. Furtado, and A. Korolev.  Mixed-phase clouds in a turbulent environment. Part 2: Analytic treatment.  Quarterly Journal of the Royal Meteorological Society, 140(680):870–880, 2014.
  * [8] C. Gardiner.  Stochastic Methods: A Handbook for the Natural and Social Sciences.  Springer-Verlag Berlin, Heildelberg, 2009.
  * [9] A. A. Hill, P. R. Field, K. Furtado, A. Korolev, and B. J. Shipway.  Mixed-phase clouds in a turbulent environment. Part 1: Large-eddy simulation experiments.  Quarterly Journal of the Royal Meteorological Society, 140(680):855–869, 2014.
  * [10] A. P. Khain and M. Pinsky.  Physical Processes in Clouds and Cloud Modeling.  Cambridge University Press, 2018.
  * [11] V. I. Khvorostyanov and J. A. Curry.  Toward the theory of stochastic condensation in clouds. part i: A general kinetic equation.  Journal of the Atmospheric Sciences, 56(23):3985 – 3996, 1999.
  * [12] A. Korolev and P. R. Field.  The effect of dynamics on mixed-phase clouds: Theoretical considerations.  Journal of the Atmospheric Sciences, 65(1):66–86, 2008.
  * [13] A. V. Korolev and I. P. Mazin.  Supersaturation of water vapor in clouds.  Journal of the Atmospheric Sciences, 60(24):2957–2974, 2003.
  * [14] G. A. Pavliotis.  Stochastic Processes and Applications, volume 60.  Springer, New York, 2014.
  * [15] T. V. Prabha, A. Khain, R. S. Maheshkumar, G. Pandithurai, J. R. Kulkarni, M. Konwar, and B. N. Goswami.  Microphysics of premonsoon and monsoon clouds as seen from in situ measurements during the cloud aerosol interaction and precipitation enhancement experiment (caipeex).  Journal of the Atmospheric Sciences, 68(9):1882 – 1901, 2011.
  * [16] P. Prabhakaran, A. S. M. Shawon, G. Kinney, S. Thomas, W. Cantrell, and R. A. Shaw.  The role of turbulent fluctuations in aerosol activation and cloud formation.  Proceedings of the National Academy of Sciences, 117(29):16831–16838, 2020.
  * [17] H. Risken.  The Fokker-Planck Equation.  Springer, second edition, 1989.
  * [18] H. C. Rodean.  Stochastic Lagrangian models of turbulent diffusion.  Meteorological Monographs ; 48. American Meteorological Society, Boston, Massachusetts, 1st ed. 1996. edition, 1996.
  * [19] R. R. Rogers and M. K. Yau.  A Short Course in Coud Physics.  Pergamon Press, third edition, 1989.
  * [20] G. Sardina, F. Picano, L. Brandt, and R. Caballero.  Continuous growth of droplet size variance due to condensation in turbulent clouds.  Physical review letters, 115(18):184501, 2015.
  * [21] R. A. Shaw.  Particle-turbulence interactions in atmospheric clouds.  Annual Review of Fluid Mechanics, 35(1):183–227, 2003.
  * [22] H. Siebert and R. A. Shaw.  Supersaturation fluctuations during the early stage of cumulus formation.  Journal of the Atmospheric Sciences, 74(4):975 – 988, 2017.
  * [23] C. Siewert, J. Bec, and G. Krstulovic.  Statistical steady state in turbulent droplet condensation.  Journal of Fluid Mechanics, 810:254–280, 2017.
  * [24] P. Squires.  The growth of cloud drops by condensation.  Australian Journal of Chemistry, 1952.
  * [25] G. E. Uhlenbeck and L. S. Ornstein.  On the theory of the brownian motion.  Phys. Rev., 36:823–841, Sep 1930.
  * [26] M. D. Zelinka, T. A. Myers, D. T. McCoy, S. Po-Chedley, P. M. Caldwell, P. Ceppi, S. A. Klein, and K. E. Taylor.  Causes of higher climate sensitivity in cmip6 models.  Geophysical Research Letters, 47(1):e2019GL085782, 2020.