%\documentclass[aps,prl,twocolumn,groupedaddress]{revtex4-1}
\documentclass[aps,pre,reprint,superscriptaddress,floatfix]{revtex4}
\usepackage{graphicx,natbib,bm,url,color}
\input macros
\graphicspath{{./fig/}{./png/}}
\newcommand{\obs}[1]{\textcolor{red}{#1}}
\newcommand{\obsb}[1]{\textcolor{blue}{#1}}

\begin{document}

\title{Cloud droplets growth by supersaturation fluctuation}
%\author{}
%\affiliation{}
\date{\today,~ $ $Revision: 1.84 $ $}
\begin{abstract}
\end{abstract}
\maketitle

\section{Introduction}

Recent studies show that supersaturation fluctuation may play
an important role in broadening the size distribution of the cloud
droplets. The laboratory study by \citet{Chandrakar16}
suggested that supersaturation fluctuations in the low aerosol/slow
microphysics limit are likely of leading importance for precipitation
formation. Field observations support the notion that both aerosol
activation and droplet growth take place in the presence of a broad
distribution of supersaturation \cite{Siebert17}.
The condensational growth due to the supersaturation fluctuation is likely
more sensitive to the large-scale properties of turbulence \cite{Gotzfried17}.
However, the coagulational growth depends on the small scales of turbulence.
When the two processes are combined, it is challenging to have large scale
separations in the simulation.

\section{Numerical model and initial setup}

We combined the Eulerian description of the turbulent velocity, temperature,
and water vapor content field with the Lagrangian description of the ensemble of
cloud droplets.

\subsection{Equations of motion for flow fields}
The turbulent velocity field $\bm{u}(\xx, t)$, the temperature field $T(\xx, t)$ and the vapor
content field $q_v(\xx, t)$ are the necessary fields that describe the moisture turbulent convection.
The Eulerian equations in the Boussinesq approximation are given by,
\begin{equation}
\label{eq:continuity}
{\partial\rho\over\partial t}+{\bm{\nabla}}\cdot(\rho\mbox{\boldmath $u$})=0 ,
\end{equation}
\begin{equation}
{\partial\mbox{\boldmath $u$}\over\partial t}+\mbox{\boldmath $u$}\cdot{\bm{\nabla}}\mbox{\boldmath $u$}={\bm f}
-\rho^{-1}{\bm{\nabla}} p
+\rho^{-1} {\bm \nabla} \cdot (2 \nu \rho \mbox{\boldmath ${\sf S}$})+B\mbox{\boldmath $e$}_z ,
\label{turb}
\end{equation}
\EQ
\label{eq:tt}
{\partial T\over\partial t}+\uu\cdot\nab T
=\kappa\nabla^2 T+\frac{L}{c_p}C_d ,  
\label{temperature}
\EN
\EQ
{\partial q_v\over\partial t}+\uu\cdot\nab q_v=D\nabla^2 q_v-C_d, 
\label{mixingRatio}
\EN
where ${\bm f}$ is a random forcing function \citep{BD02},
$\nu$ is the kinematic viscosity of the air flow,
${\sf S}_{ij}={\textstyle\frac{1}{2}}(\partial_j u_i+\partial_i u_j)
-{\textstyle{1\over3}}\delta_{ij}{\bm{\nabla}}\cdot\mbox{\boldmath $u$}$ is
the traceless rate-of-strain tensor, $p$ is the gas pressure, $\rho$ is the gas density,
$c_p$ is the specific heat at constant pressure, and $L$ is
the latent heat. $\kappa$ is the thermal diffusivity coefficient of the air and $D$
is the water vapor diffusivity coefficient.
We assume that the gas is isothermal with constant sound speed $c_{\rm s}$.
The buoyancy field depends on the temperature field $T$, vapor mixing ratio $q_v$,
and the liquid mixing ratio $q_l$, 
\EQ
%B(\mbox{\boldmath $x$},t)=g[\frac{T-T_a}{T_a}+\epsilon(q_v-q_{va})-q_l],
B(\xx,t)=g[\frac{T^{\prime}}{T}+\epsilon q_v^{\prime}-q_l],
\EN
where $\epsilon=R_v/R_d-1\approx0.608$. Here, $R_v$ and $R_d$ are the vapor
and dry air gas constant respectively.
The temperature fluctuations are given by
\EQ
T^{\prime}\left(\xx,t\right)=T\left(\xx,t\right)-\left<T(t)\right>,
\EN
and the vapour mixing ratio fluctuations by
\EQ
q_v^{\prime}\left(\xx,t\right)=q_v\left(\xx,t\right)-\left<q_v(t)\right>,
\EN
with $\left<\cdot\right>$ is a volume average \cite{Gotzfried17}.
To simplify the calculation, we adopt the same method in \cite{Kumar14},
where the mean environmental temperature and vapor mixing ratio are
constant. This assumption is plausible in the circumstance that we
do not consider the entrainment, i.e., there is only mass and energy
transfer between liquid water and water vapor.
The condensation rate $C_d$ is given by \cite{Kumar14, saito2017turbulence} 
\EQ
C_d\left(\xx,t\right)=
\frac{1}{m_a}\frac{dm_l\left(\xx,t\right)}{dt}= 
\frac{4\pi\rho_l G}{\rho_a (\Delta x)^3}\sum_{j=1}^{N_{\triangle}}  
s\left(\xx,t\right)r\left(t\right). 
\label{condensationRate}
\EN
where $m_a$ and $m_l$ are the mass of 
air and liquid water, respectively,
$\rho_l$ and $\rho_a$ are the water density 
and the reference mass density of dry 
air, respectively, $N_{\triangle}$ is the total 
number of droplets in the grid cell $(\Delta x)^3$. 
To conserve the mass of liquid water, 
we define the mixing ratio of total water 
as $q_{tot}=q_l+q_v+q_a$, where $q_l$ is the mixing 
ratio of liquid water defined as \cite{saito2017turbulence},
\EQ
q_l\left(\xx,t\right)= 
\frac{4\pi\rho_l}{3\rho_a (\Delta x)^3}\sum_{j=1}^{N_{\triangle}}  
r\left(t\right)^3. 
\EN

The supersaturation rate $s$ is defined as 
\EQ
s\left(\xx,t\right)=\frac{q_v\left(\xx,t\right)}{q_{vs}\left(T\right)}-1, 
\label{supersaturation}
\EN
where $q_{vs}\left(T\right)$ is the 
supersaturated vapor pressure at 
temperature $T$ and can be determined via the
Clausius-Clapeyron equation. It is given by
\EQ
q_{vs}\left(T\right)=\frac{e_s\left(T\right)}{R_v \rho_0 T}. 
\label{eq: qvs}
\EN
$R_v$ is the specific gas constant for water vapour $R_v$. The saturation
pressure $e_s$ is given by
\EQ
\label{eq:e_s}
e_s(T)=c_1 \exp(-c_2/T),
\EN
where $c_1$ and $c_2$
are constants. We refer to Table 2 in \cite{Gotzfried17} for the constants.

\subsection{Lagrangian model for cloud droplets}
The superparticle \cite{li16} is treated as a
Lagrangian point-particle, where one solves for the particle position
$\mbox{\boldmath $x$}_i$ via
\begin{equation}
\frac{d\mbox{\boldmath $x$}_i}{dt}=\mbox{\boldmath $V$}_i
\end{equation}
and the velocity via
\begin{equation}
\frac{d\mbox{\boldmath $V$}_i}{dt}=\frac{1}{\tau_i}(\mbox{\boldmath $u$}-\mbox{\boldmath $V$}_i)+\bm{g}
\end{equation}
in the usual way.
Here, $\mbox{\boldmath $u$}$ is the fluid velocity at the position of the superparticle,
$\bm g$ is the gravitational acceleration, 
$\tau_i$ is the particle inertial response or stopping time
of a particle in superparticle $i$ and is given by
\begin{equation}
\tau_i=\frac{2\rho_d r_i^2}{9\rho\nu_i^{\rm eff}},
\label{stopping_time}
\end{equation}
where $r_i$ is the radius of particles in superparticle $i$,
$\rho_{\rm d}$ is the particle solid material density 
and the effective viscosity is given by 
\begin{equation}
\nu_i^{\rm eff}=\nu\,(1+0.15\,{\rm Re}_i^{0.687}),
\end{equation}
where $\nu$ is the ordinary (microphysical) fluid viscosity, and
${\rm Re}_i=2r_i|\mbox{\boldmath $u$}-\mbox{\boldmath $V$}_i|/\nu$ is the particle Reynolds number,
which provides a correction factor to the particle stopping time.

The growth of the particle radius $r_i$ by condensation
is governed by \cite{2011_lamb,2015_Sardina}
\begin{equation}
	{{\rm d} r_i\over{\rm d} t}={Gs\left(\xx_i,t\right)\over r_i},
\label{cond_eq}
\end{equation}
where $s$ is the supersaturation given by \Eq{supersat}
and $G$ is the condensation parameter
(having units of $\,{\rm m}^2\,{\rm s}^{-1}$). 

\section{Result}

We conduct high-resolution simulations for different cases.
Table~\ref{Swarm_Rey} shows the initial time scales.
The initial configuration is unsaturated with $s=-0.0168$.

\begin{table}[t!]
\caption{Summary of the simulations.}
\centering
\setlength{\tabcolsep}{3pt}
\begin{tabular}{lr}
Quantity ($\rm{s}$) & DNS \\
\hline
$\tau_{\eta}$ & 0.16 \\ 
$\tau_{\rm{p}}$ & 0.002 \\ 
$\tau_{\rm{ph}}=(4\pi D^{\prime}n_0r_{\rm ini})^{-1}$ & 9.448 \\ 
$\tau_{\rm{eddy}}$ & 0.246 \\ 
$\tau_{\rm{ev}}=r_{\rm ini}^2/(2Gs_0)$ & 25.433 \\ 
\hline
\multicolumn{2}{p{0.48\textwidth}}{
Here, $n_0=1.25\times10^8\,\rm{m}^{-3}$ and $r_{\rm ini}=10\,\mu$m.
$D^{\prime}=6.74\times10^{-6}\,\rm{m}^2\rm{s}^{-1}$ is the modified vapor diffusivity,
and $G=1.17\times10^{-10}\,\rm{m}^2\rm{s}^{-1}$ (we adopt values from Table 2 of \cite{Gotzfried17}). 
}
%3D/python/initial_time.py
\end{tabular}
\label{Swarm_Rey}
\end{table}

\Fig{timeSeries_SW128condensP_coag0_grav0_buoyancy0_n1e8}
shows the time series of the field quantities.
The fluctuations of $q_v$ and $T$ does not change much
compared with the mean value. The condensation rate
$C_d$ and the supesaturation ratio $s$ fluctuates dramatically. 

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{timeSeries_SW128condensP_coag0_grav0_buoyancy0_n1e8}
\end{center}\caption[]{Time series of the field quantities.
$\nu=1.5\times10^{-5}\m^2\s^{-1}$, $\kappa=2.55\times10^{-5}\m^2\s^{-1}$ so that $Pr=0.4$,
and $D=2.52\times10^{-5}\m^2\s^{-1}$ so that the Schmit number $Sc=\nu/D=0.4$.
$L=0.125\m$. $q_v(\xx, t=0)=0.0163\,kg \cdot kg^{-1}$,
$\nabla q_0=(0.01, 0.01, 0.01)$, $T_0(\xx, t=0)=293.25\, K$ and $\nabla T_0=(100, 100, 100)$.
The environmental temperature is $T_e=293.00\,K$ and
the environmental water vapor mixing ratio is $q_{v,e}(\xx, t=0)=0.0163\,kg \cdot kg^{-1}$. 
A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted. $n_0=1.25\times10^8\m^{-3}$.
$N_{\rm s}=0.058N_{\rm grid}$. 
$G=1.17\times 10^{-10}\m^2\s^{-1}$.
}
\label{timeSeries_SW128condensP_coag0_grav0_buoyancy0_n1e8}
\end{figure*}

The mixing of the clear air and cloud air can be characterized
by the Damk\"{o}hler number, the ratio of a fluid time scale
to a characteristic thermodynamic time scale associated with
the evaporation process
\EQ
Da=\frac{T_L}{\tau_{\rm phase}}.
\label{Damkohler}
\EN
$Da\ll1$ characterizes the homogeneous mixing and $Da\gg1$ characterizes
the inhomogeneous mixing. $\tau_{\rm phase}$ is defined in the appendix.
\Fig{DamKohler_SW128condensP_coag0_grav0_buoyancy0_n1e8} shows the time
evolution of $Da$. $Da\ll1$ in our simulation. Therefore, the homogeneous
mixing dominates the condensation.

\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{DamKohler_SW128condensP_coag0_grav0_buoyancy0_n1e8}
\end{center}\caption[]{Time evolution of $Da$.
Same simulation as \Fig{timeSeries_SW128condensP_coag0_grav0_buoyancy0_n1e8}.} 
\label{DamKohler_SW128condensP_coag0_grav0_buoyancy0_n1e8}
\end{figure}

Comparison among different cases is shown in \Fig{f_grav0Re50_cond_coa_comp}.
Condensation due to supersaturation fluctuations does not contribute
much to the broadening of the size distribution. This is because that
$n_0=10^{10}\m^3$ is about 100 times larger than $n_{\rm ref}=10^8\m^3$.
When $n_0=10^8\m^3$, the broadening becomes outstanding as shown in
\Fig{f_grav0Re50_cond_coa_comp_n1e8}. In the present study, the updraft
velocity of the simulation box is zero. Therefore, the condensational
growth of cloud droplets is dominated by the supersaturation fluctuations.
This can be used to model the stratocumulus cloud, where the updraft velocity
of the parcel is almost zero \cite{korolev1995influence}.


\begin{figure}[h!]\begin{center}
%\includegraphics[width=0.49\textwidth]{m_grav0Re50_cond_coa_comp}
\includegraphics[width=0.49\textwidth]{f_grav0Re50_cond_coa_comp}
\end{center}\caption[]{Comparison among different cases.
$\nu=10^{-5}\m^2\s^{-1}$, $\kappa=2.55\times10^{-5}\m^2\s^{-1}$ so that $Pr=0.4$,
and $D=2.52\times10^{-5}\m^2\s^{-1}$ so that the Schmit number $Sc=\nu/D=0.4$.
$L=0.125\m$. $q_v(\xx, t=0)=0.0163\,kg \cdot kg^{-1}$,
$\nabla q_0=(0.01, 0.01, 0.01)$, $T_0(\xx, t=0)=293.25\, K$ and $\nabla T_0=(100, 100, 100)$. 
A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted. $n_0=10^{10}\m^{-3}$.
$N_{\rm s}=N_{\rm grid}/2$. 
$G=1.17\times10^{-10}\m^2\s^{-1}$.
}
\label{f_grav0Re50_cond_coa_comp}
\end{figure}

\begin{figure}[h!]\begin{center}
%\includegraphics[width=0.49\textwidth]{m_grav0Re50_cond_coa_comp}
\includegraphics[width=0.49\textwidth]{f_grav0Re50_cond_coa_comp_n1e8}
\end{center}\caption[]{Same as \Fig{f_grav0Re50_cond_coa_comp} but with $n_0=10^8\m^{-3}$.} 
\label{f_grav0Re50_cond_coa_comp_n1e8}
\end{figure}

%For other simulations, we may adopt unsaturated initial configuration.
The initial configuration of the mixing ratio and temperature
are given as the same as the one in \cite{Gotzfried17}, where the values are
matched to the ones obtained in the CARRIBA experiments \cite{katzwinkel2014measurements}.
This configuration means the water vapor is saturated initially.

Following our previous study \cite{li2017eulerian} of the collisional
growth, we investigate the $Re_{\lambda}$ and $\bar{\epsilon}$ dependencies
of the condensational growth as shown in \Fig{f_gravRe50_cond_e_comp2}.
The details of the simulations are shown in Table~\ref{Swarm_Rey}
\begin{table*}[t!]
\caption{Summary of the simulations.}
\centering
\setlength{\tabcolsep}{3pt}
\begin{tabular}{lcccccccccr}
	Run &  $f_0$ & L (\rm{m}) & $N_{\rm grid}$ & $N_{\rm s}$& $u_{\rm rms}$ ($\rm{m} \, \rm{s}^{-1}$) & $\mbox{\rm Re}_{\lambda}$ & $\bar{\epsilon}$ ($\rm{m}^2\rm{s}^{-3}$)& $\eta\cdot 10^{-4}$ ($\rm{m}$) & $\tau_{\eta}$ ($\rm{s}$) & $\tau_{\rm eddy}$ ($\rm{s}$) \\
\hline
A & $0.02$  & 0.125 & $128^3$ & $N_{\rm s, 128}$ &0.16 &45  & 0.037 & 5.5  & 0.02 & 0.25 \\ %3D/SW128condensP_coag0_grav_buoyancy_n2p5e8
B & $0.02$  & 0.25  & $256^3$ & $2^3N_{\rm s, 128}$&0.22 &79  & 0.039  & 5.4  & 0.02 & 0.37 \\ %3D/SW256condensP_coag0_grav_buoyancy_n2p5e8_L0p25
C & $0.007$ & 0.164 & $128^3$ & $N_{\rm s, 128}$&0.11 &42 & 0.008 & 8.2  & 0.017 & 0.50 \\ %3D/SW128condensP_coag0_grav_buoyancy_n2p5e8_f7em3_L0p164
D & $0.002$ & 0.222 & $128^3$ & $N_{\rm s, 128}$&0.13 &61 & 0.008 & 8.1  & 0.044 & 0.56 \\ %3D/SW128condensP_coag0_grav_buoyancy_n2p5e8_f2em3_L0p222
$D^{\ast}$ & $0.007$ & 0.222 & $128^3$ & $N_{\rm s, 128}$&0.14 &61 & 0.01 & 7.6  & 0.038 & 0.52 \\ %3D/SW128condensP_coag0_grav_buoyancy_n2p5e8_f7em3_L0p222
E & $0.02$  & 0.5  & $512^3$ & $2^6N_{\rm s, 128}$&0.32 &142  & 0.059  & 4.9  & 0.016 & 0.50 \\ %3D/SW512condensP_coag0_grav_buoyancy_n2p5e8_L0p5
F & $0.05$  & 0.1 & $128^3$ & $N_{\rm s, 128}$ &0.28 &55  & 0.231 & 3.5  & 0.008 & 0.11 \\ %3D/SW128condensP_coag0_grav_buoyancy_n2p5e8_f0p05_L0p1
F1 & $0.05$  & 0.08 & $128^3$ & $N_{\rm s, 128}$ &0.26 &45  & 0.237 & 3.5  & 0.008 & 0.10 \\ %3D/SW128condensP_coag0_grav_buoyancy_n2p5e8_f0p05_L0p08
G & $0.014$  & 0.15 & $128^3$ & $N_{\rm s, 128}$ &0.14 &46  & 0.019 & 6.5  & 0.028 & 0.35 \\ %3D/SW128condensP_coag0_grav_buoyancy_n2p5e8_f0p014_L0p15
\hline
\multicolumn{11}{p{0.85\textwidth}}{
Here, $f_0$ is the amplitude of the random forcing (see text) and $L$ is the domain size.
$\nu=1.5\times10^{-5}\m^2\s^{-1}$. The liquid water content is fixed, thus $n_0=2.5\times10^8\m^{-3}$
is same for all the simulations. To avoid artificial acceleration of the collision rate, we
keep $n_{\rm p/s}$ the same since the collision rate is proportional to $n_{\rm p/s}$.
When $L=0.125$, $N_{\rm s}$ is determined such that $N_{\rm p}/N_{\rm s}=2$.
This results in $N_{\rm s}/N_{\rm grid}\approx0.1$, which is still within the convergence
range $N_{\rm s}/N_{\rm grid}\approx0.05$. Another important reason to keep $N_{\rm s}/N_{\rm grid}\approx0.1$
is that in a cloud system, there are about 0.1 droplet per cubic Kolmogorov scale (assume the length of
the grid cell is the same as $\eta$).
We denote $N_{\rm s, 128}=244140$. When varying $N_{\rm grid}$, $N_{\rm s}=2^{\beta} N_{\rm s, 128}$,
where $\beta$ is an integer.
}
\end{tabular}
\label{Swarm_Rey}
\end{table*}

\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{f_gravRe50_cond_e_comp2}
\end{center}\caption[]{Comparison of the size distribution for
different $Re_{\lambda}$ and $\bar{\epsilon}$. Buoyancy force and
gravity are all turned on.
$\nu=1.5\times10^{-5}\m^2\s^{-1}$, $\kappa=\nu$ so that $Pr=1.$,
and $D=2.55\times10^{-5}\m^2\s^{-1}$ so that the Schmit number $Sc=\nu/D=0.6$.
$q_v(\xx, t=0)=0.0157\,kg \cdot kg^{-1}$,
$\nabla q_0=(0.01, 0.01, 0.01)$, $T_0(\xx, t=0)=292\, K$ and $\nabla T_0=(100, 100, 100)$. 
$q_{v,en}=0.01\,kg \cdot kg^{-1}$ and $T_{v,en}=293\,K$. 
A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted. $n_0=2.5\times10^8\m^{-3}$.
$N_{\rm s}/N_{\rm grid}\approx0.1$. 
$G=1.17\times10^{-10}\m^2\s^{-1}$.
We use the same initial configurations for thermodynamics as the ones
in \cite{Gotzfried17}. The initial mean supersaturation is $\bar{s}=0.02$,
meaning the initial environment is saturated.
See details of the simulations in Table~\ref{Swarm_Rey}.
} 
\label{f_gravRe50_cond_e_comp2}
\end{figure}

We further investigate $\tau{\rm phase}$ in \Fig{DamKohler_gravRe50_cond_e_comp}.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{DamKohler_gravRe50_cond_e_comp}
\end{center}\caption[]{$Da$ and $\tau_{\rm phase}$.
Same simulations as in \Fig{f_gravRe50_cond_e_comp2}.
} 
\label{DamKohler_gravRe50_cond_e_com}
\end{figure*}

The corresponding time series is shown in
\Fig{timeSeries_gravRe50_cond_e_comp}.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{timeSeries_gravRe50_cond_e_comp}
\end{center}\caption[]{Time series.
Same simulations as in \Fig{f_gravRe50_cond_e_comp2}.
} 
\label{timeSeries_gravRe50_cond_e_comp}
\end{figure*}

We also check the spectra of the number density of particles
in \Fig{Ec_gravRe50_cond_e_comp}. This aims to investigate the fluctuations of
$\tau{\rm phase}$.

\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{Ec_gravRe50_cond_e_comp}
\end{center}\caption[]{Spectra of the number density of particles.
Same simulations as in \Fig{f_gravRe50_cond_e_comp2}.
} 
\label{Ec_gravRe50_cond_e_comp}
\end{figure}

Combined condensational and collision growth is shown in
\Fig{f_evolution_cond_coa_comp}.

\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{f_evolution_cond_coa_comp}
\end{center}\caption[]{Combined condensation and collisional growth.
See details of the simulation in Run B in the table of the paper.} 
\label{f_evolution_cond_coa_comp}
\end{figure}

\appendix

\Tab{constants} shows the list of constants for the thermodynamics used in the present paper.
\label{tab:constant}

\begin{table}[t!]
\caption{List of constants for the thermodynamics.}
\centering
\setlength{\tabcolsep}{3pt}
\begin{tabular}{lr}
Quantity & Value \\
\hline
$\nu$ ($\m^2\s^{-1}$) & $1.5\times10^{-5}$ \\ 
$\kappa$ ($\m^2\s^{-1}$) & $1.5\times10^{-5}$ \\ 
$D$ ($\m^2\s^{-1}$) & $2.55\times10^{-5}$ \\ 
$G$ ($\m^2\s^{-1}$) & $1.17\times10^{-10}$ \\ 
$c_1$ (Pa) & $2.53\times10^{11}$ \\ 
$c_2$ ($K$) & 5420 \\ 
$L$ ($J \cdot kg^{-1}$) & $2.5\times10^6$ \\ 
$c_{\rm p}$ ($J \cdot kg^{-1} K^{-1}$) & 1005.0 \\ 
$R_v$ ($J \cdot kg^{-1} K^{-1}$) & 461.5 \\ 
$\alpha$ & 0.608 \\
$Pr=\nu/\kappa$ & 1 \\
$Sc=\nu/D$ & 0.6 \\
$q_v(\xx, t=0)$ ($\,kg \cdot kg^{-1}$) & 0.0157 \\
$\nabla q_0$ ($\,kg \cdot kg^{-1}$) & ($10^{-4}$, $10^{-4}$, $10^{-4}$)\\
$q_{v,\rm{en}}\,(kg \cdot kg^{-1})$ & 0.01\\
$T_0(\xx, t=0)\,(K)$ & 292 \\
$\nabla T_0$ ($K$) & (1, 1, 1) \\ 
$T_{v,\rm{en}}\,(K)$ & 293 \\ 
\hline
\multicolumn{2}{p{0.48\textwidth}}{
This initial configuration gives a mean supersaturation of $\bar{s(\xx,t=0)}=0.02$,
meaning the initial environment is saturated.
%XY: notes. Will be removed when submitting the manuscript.
%We adopt values from Table 2 of \cite{Gotzfried17}. 
%$\Delta_z qv=10^{-5}\,kg\cdot kg^{-1}$ is adopted from the observation \cite{osman2018characteristics}.
%$\Delta_z T=0.01\,K$ is employed from the observation \cite{siebert2017supersaturation}.
%The rms values from the DNS of \cite{2015_Sardina} are $s_{\rm rms}=0.006$ (C1) and $s_{\rm rms}=0.0075$ (C2)
%when $Re_{\lambda}=150$.
%XY.
}
%3D/python/initial_time.py
\end{tabular}
\label{constants}
\end{table}
\section{Generalized Twomey model of supersaturation}
The spatial and temporal evolution of the supersaturation
within the clouds can can be described as \cite{2011_lamb}
(P417, Ch10; P248, Ch6)
\begin{equation}
\frac{Ds}{Dt}=A_1w-A_2\frac{dw_L}{dt}+\kappa\Delta s
\end{equation}
where $s=s(\xx,t)$ is the supersaturation, 
$w_L$ the liquid water mass content, 
$A_1=5\times10^{-4}\m^{-1}$, 
$A_2=350\kg^{-1}\m^3$, 
$A_3=5\times10^{-11}\m^2\s^{-1}$ \cite{2009_Lanotte}, 
the thermal parameter. It is worth noting that here $s$
is a mean field quantity. When it is applied to \Eq{cond_eq},
mapping from Eulerian field to Lagrangian field is required.
$\kappa \approx 10^{-5}\rm{m}^2\rm{s}^{-1}$ is
the thermal diffusion coefficient of water vapor, $w$ is the
vertical component of the velocity of the air flow.
For a quasi-steady state \cite{2011_lamb} (Eq(10.13), P424, Ch10),
\begin{equation}
\frac{dw_L}{dt} = \int_{0}^{\infty} \frac{dm_d}{dt}f_d dr_d
= 4 \pi \rho_w A_2 A_3 s \int_{0}^{\infty} r_d f_d dr_d
\end{equation}

Define relaxation time \cite{2011_lamb} (Eq(10.17), P425, Ch10):  
\begin{equation}
\tau_s^{-1}(\xx,t) = 4 \pi \rho_w A_2 A_3 \int_{0}^{\infty} r_d f_d dr_d,
\end{equation}
we adopt the nomenclature in the paper \cite{li16}, the relaxation
time can be expressed as $\tau_s^{-1} = 4 \pi \rho_w A_2 A_3\bar{r}M_0$
or alternatively, as $\tau_s^{-1} = 4 \pi \rho_w A_2 A_3M_1$.
The relaxation time is inversely proportional to the mean radius and
mean number density, indicating that larger cloud droplets or droplets
with high number density can ``eat up'' the water vapor quickly so that
the system can restore to equilibrium state faster.
The discretized form of the relaxation time can be written as
\begin{equation}
	\tau_s^{-1}=4 \pi \rho_w A_2 A_3 \sum_{k=1}^{r_{\rm max}} f_k r_k\delta r_k. 
\label{relaxation}
\end{equation}
The phase transition time is determined by the first moment
of the size distribution. Vice versa, the size distribution
is also governed by the phase transition time. \cite{desai2018influence}
suggested that the variations of the mean radius and the mean
number density contribute significantly to the increase rate
of the size distribution width, not to the mean growth. It is
therefore important to study the source of the variability,
for example, particle clustering due to turbulence and gravity.
We will compare the power spectra of the number density $E_n$
for different $Re_{\lambda}$ and $\bar{\epsilon}$.

Thus, the supersaturation field $s$ is transported 
by the fluid by \cite{2009_Lanotte,2015_Sardina}: 
\begin{equation}
\frac{Ds}{Dt}=A_1w+\kappa\Delta s - \frac{s}{\tau_s},
\label{supersat}
\end{equation}
which can also be derived from first principles, vapor
mass and energy conservation, and thermodynamics relations;
see the supplemental material of the paper \cite{2015_Sardina}. 
Note that the definition of supersaturation \cite{2011_lamb,2015_Sardina}
(Eq(3.73), P115, Ch3 and P321, Ch8) is
\EQ
s=\frac{\Delta e_v}{e_{vs}}=\frac{e_v-e_{vs}}{e_{vs}}=\frac{e_v}{e_{vs}}-1, 
\EN
which is different from ``saturation ratio'' $S=s+1$.
The supersaturation can also be expressed as a function
of vapor mixing ratio by adopting ideal gas law,
\EQ
s=\frac{e_v}{e_{vs}}-1=\frac{\rho_vR_vT}{\rho_{vs}R_vT}-1=\frac{\rho_v/\rho_d}{\rho_{vs}/\rho_d}-1=\frac{q_v}{q_{vs}}-1,
\EN
where $q_v$ is the vapor mixing ratio defined as $q_v=\rho_v/\rho_d$
and $q_{vs}$ is the saturated vapor mixing ratio.
When solving \Eq{supersat}, the vertical velocity $w$ should be the velocity
of cloud droplets. Some studies \cite{2015_Sardina} obtain $w$ by interpolating
from the Eulerian grid to Lagrangian points. $\tau$ is evaluated at first at
the droplet position and then redistributed to the nearest eight grid vertices \cite{Kumar2013,2015_Sardina}.

%\Fig{supersat} shows simple tests about the passive
%scalar equation.
%
%\begin{figure}[h!]\begin{center}
%\includegraphics[width=0.49\textwidth]{supersat}
%\end{center}\caption[]{To test the passive scalar equation
%of supersaturation in 1-D turbulence case. $A_3=5\times10^{-15}\m^{-1}$
%so that the code does not blow up.
%%particles_coag/1D/SW16condensP_coag0_c
%}
%\label{supersat}
%\end{figure}

\section{Size distribution}
To calculate $\tau_s$ in a simple manner, we need to understand
the size distribution of cloud droplets. The number density distribution
function (size distribution in meteorological context) is initially
from statistical physics. In the phase space (kinetic theory),
it is a function of
position, momentum, size, and time: $f(\xx,\vv_{\rm phase},r,t)$.
In the present study, we simplify it to a hydrodynamic quantity
as $f(\xx,r,t)$. It is defined as
\EQ
f(\xx,r,t)=\lim_{V\to0} \frac{dN}{dD},
\EN
where $V$ is the volume at position $\xx$, 
$dN$ is the number of particles with size
$D$ between size range $dD$ in the volume $V$. When calculating
$f(\xx,r,t)$ at each time step in DNS , we need to choose $V$ properly
to represent $f$ at position $\xx$. 
\section{Implementation of the supersaturation in the Pencil-Code}
Since the supersaturation is transported by background air flow
and interacted with cloud droplets, we need to map the source term
$A_1w$ and the sink term $s/\tau_s$ to the same grid of tracer
particles. Thus we re-write \Eq{supersat} as  
\begin{equation}
	\frac{Ds}{Dt}=A_1w+\kappa\Delta s
	- \frac{s}{\tau_s}\delta(\xx_p-\XX),
\label{supersat_numerical}
\end{equation}
In the Pencil-Code, at fixed position $\xx$, 
we have $f(\xx,r,t)$ and $r$. Therefore, the only
thing we need to do is to map the Lagrangian field
to Euler field. With a fixed Euler grid cell, $\tau_s$
can be calculated by \Eq{relaxation}.

%First, we test the implementation of the diffusion equation
%(\Eq{supersat} without the updraft and relaxation terms).
%As shown in \Fig{ssatrms_SW16condensP_coag0c}, 
%
%\begin{figure}[h!]\begin{center}
%\includegraphics[width=0.49\textwidth]{ssatrms_SW16condensP_coag0c}
%\end{center}\caption[]{Test without updraft and relaxation terms.
%The diffusion coefficient is set to $\kappa=0.1\,\rm{m}^2\rm{s}^{-1}$.
%}
%\label{ssatrms_SW16condensP_coag0c}
%\end{figure}

\section{Spikes of $dt$}

\Fig{condens_dt} shows that the wrong spikes.
\Fig{condens_dt2} shows the correct result with new implementation
of the time step contribution. However, for high resolution simulations,
the spikes appear again. Therefore, we debug the code using high resolution
simulations. \Fig{condens_dt3} demonstrates the results with
$\nu=1.5\times10^{-5}\m^2\s^{-3}$.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_dt}
\end{center}\caption[]{Spikes of $dt$.
$\nu=5\times10^{-4}\m^2\s^{-1}$, $\kappa=\nu$ so that
the Prandtl number is $Pr=\nu/\kappa=1$, and $D=\nu$ so that the Schmit number $Sc=\nu/D=1$.
$N_{\rm grid}=16^3$, $N_{\rm s}=N_{\rm grid}/2$.
$L=0.1\m$.
A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted. $n_0=10^{10}\m^{-3}$.
The scheme for the time step contribution is from Anders.
}
\label{condens_dt}
\end{figure*}

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_dt2}
\end{center}\caption[]{Same simulations as in \Fig{condens_dt} but the
time step contribution is different, which is now based on $\tau_{\rm cond}$
and $\tau_{\rm evap}$.
}
\label{condens_dt2}
\end{figure*}

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_dt3}
\end{center}\caption[]{Same simulations as in \Fig{condens_dt2} but with $\nu=1.5\times10^{-5}\m^2\s^{-1}$.
$\kappa=5\times10^{-4}\m^2\s^{-1}$ and $D=\nu$. $n_0=10^8\m^{-3}$.
The simulation with $G=1.17\times10^{-10}\m^2\s^{-1}$ crashes as expected.
}
\label{condens_dt3}
\end{figure*}

\Fig{condens_dt_N256} shows the simulation with or without particles.
However, when particle modules are turned on, the time steps exhibit
spikes with or without condensation. To investigate the spikes, we conducted
simulations with different resolution. As shown in \Fig{condens_dt_resolution},
higher resolution reveals more spikes.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_dt_N256}
\end{center}\caption[]{Spikes of $dt$ for high resolution simulation.
$\nu=1.5\times10^{-5}\m^2\s^{-1}$, $\kappa=\nu$ so that
the Prandtl number is $Pr=\nu/\kappa=1$ (in some cases, $\kappa=10^{-4}\m^2\s^{-1}$,
this leads to different $T_{\rm rms}$), and $D=\nu$ so that the Schmit number $Sc=\nu/D=1$.
$N_{\rm grid}=256^3$, $L=0.125\m$. $q_v(\xx, t=0)=0.0163\,kg \cdot kg^{-1}$,
$\nabla q_0=(0.01, 0.01, 0.01)$, $T_0(\xx, t=0)=293.25\, K$ and $\nabla T_0=(100, 100, 100)$. 
A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted. $n_0=10^{10}\m^{-3}$.
$N_{\rm s}=N_{\rm grid}/2$. 
The scheme for the time step contribution is based
on $\tau_{\rm evap}$ and $\tau_{\rm cond}$ or Anders's scheme. $G=0\m^2\s^{-1}$.
This demonstrates that the spikes are due to the time step contribution
of condensation in $\rm{particles\_radius.f90}$.
}
\label{condens_dt_N256}
\end{figure*}

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_dt_N256_evap_cond}
\end{center}\caption[]{Spikes of $dt$ for high resolution simulation.
The setup is the same as the one in \Fig{condens_dt_N256}. We check
different scheme of the time step contribution. Since
$\tau_{\rm phase}\ll \tau_{\eta}$, the time step contribution can be ignored.}
\label{condens_dt_N256_evap_cond}
\end{figure*}

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_dt_resolution}
\end{center}\caption[]{Spikes of $dt$ for high resolution simulation.
$\nu=1.5\times10^{-5}\m^2\s^{-1}$, $\kappa=10^{-4}\m^2\s^{-1}$,
and $D=\nu$ so that the Schmit number $Sc=\nu/D=1$.
$L=0.125\m$. $q_v(\xx, t=0)=0.0163\,kg \cdot kg^{-1}$,
$\nabla q_0=(0.01, 0.01, 0.01)$, $T_0(\xx, t=0)=293.25\, K$ and $\nabla T_0=(100, 100, 100)$. 
A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted. $n_0=10^{10}\m^{-3}$.
$N_{\rm s}=N_{\rm grid}/2$. 
The scheme for the time step contribution is based
on $\tau_{\rm evap}$ and $\tau_{\rm cond}$. $G=0\m^2\s^{-1}$.
}
\label{condens_dt_resolution}
\end{figure*}


\section{Mass conservation}

The sum of vapor and liquid water mass density in the volume $V$, $\rho_v$
and $rho_l$ has to be constant with respect to time.
As shown in \Fig{waterMassConservation}, $\rho_{\rm{total}}(t)$ is constant.
\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{waterMassConservation}
\end{center}\caption[]{Mass conservation of water.
%3D/waterMass
}
\label{waterMassConservation}
\end{figure}

\section{Reduced speed of sound technique}

Since {\sc Pencil Code} is designed to deal with weakly compressible flow,
the high temperature leads to small time step. To maintain the weak
compressibility and save the computational cost, we adopt the reduced
speed of sound technique (RSST) developed by \citet{hotta2012numerical}.
In the {sc Pencil Code}, the gas flow is assumed to be
isothermal with constant sound speed $c_s$,
so that $c_{\rm s}^2=\gamma p/\rho$, where $\gamma=c_{\rm p}/c_{\rm v}=5/3$
is the ratio between specific heats, $c_{\rm p}$ and $c_{\rm v}$, at constant
pressure and constant volume, respectively.
Using ideal gas law $p=\rho R_{\rm specific}T$, we can obtain
$c_s^2=\gamma T R_{\rm specific}$. Therefore, $c_s$ increases with
temperature. When the temperature module ($\rm{temperature\_idealgas.f90}$),
is turned on, $c_s$ is no longer constant, instead, it is determined
by the density and pressure field. If the initial temperature is too large
$c_s$ will become larger, which increases the compressibility and results
in small time step.

The RSST introduced the effective sound speed
\begin{equation}
\hat{c}_{\rm s}=\frac{1}{\xi}\sqrt{\frac{p}{\rho}}
\end{equation}
so that the continuity equation (\Eq{eq:continuity}) can be expressed as
\begin{equation}
\label{eq:continuity2}
{\partial\rho\over\partial t}+\frac{1}{\xi^2}{\bm{\nabla}}\cdot(\rho\mbox{\boldmath $u$})=0. 
\end{equation}
Using this equation, the effect speed of sound becomes $\xi$ times smaller,
but the dispersion relationship for sound waves remains otherwise unchanged
\cite{hotta2012numerical}. \citet{hotta2012numerical} found that RRST is valid
when $Ma<0.7$ and $\xi \le 20$, where the Mach number is defined as the ratio
between the rms velocity and the effective sound speed $\hat{c}_{\rm s}$.
\Fig{rrst} shows a test for different $\xi$. Apparently, the RRST does
not work correctly at the moment.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{rrst}
\end{center}\caption[]{RRST for different $\xi$.
Mach number $Ma\le 0.08$ in all the tests.
}
\label{rrst}
\end{figure*}

\section{Artificial mean temperature profile}

\subsection{Convergence of initial condition of $T$}

Since the RRST technique does not work, we employ an artificial mean temperature profile
to tackle with the compressible flow. We do not consider uplifting cooling in our simulation,
temperature fluctuations are small compared with the mean temperature.
More importantly, \citet{saito2017turbulence} found that important information is lost
if the whole temperature field is solved. Thus, they solve temperature fluctuations and
the mean temperature separately. Their temperature fluctuations were solved using \Eq{eq:tt}
and the mean temperature was solved by a parcel model.
In the present study, temperature fluctuations are solved using \Eq{eq:tt} with an initial
constant profile of $T^{\prime}_0=3.25$ K. This results in a small sound speed $c_{\rm s}$ so that
the flow is nearly incompressible.
When calculating the saturation pressure, we employ
an artificial temperature profile $\left<T\right>_0=290$ K so that
$T=T^{\prime}+\left<T\right>_0$ in \Eq{eq:e_s}. This also leads
to an initial temperature profile $T_0=293.5$ K.

As a test of the implementation, we use $\nu=5\times10^{-4}\m^2\s^{-1}$, $\kappa=\nu$ so that
the Prandtl number is $Pr=\nu/\kappa=1$, and $D=\nu$ so that the Schmit number $Sc=\nu/D=1$.
$N_{\rm grid}=16^3$, $N_{\rm s}=N_{\rm grid}/2$. A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted.
Since the flow is weakly compressible, we first check the effect of $T^{\prime}_0$ on the flow
field. When $T^{\prime}_0 \le 3.5$ K, the code crashes, which is due to the
$C_d$. To verify that the $T^{\prime}_0$ does not affect the flow field in the absence of buoyancy
force, we also conducted a simulation
without condensation. As shown in \Fig{condens_meanT_compCs}, $\urms$ and $\bar{\epsilon}$
are almost the same for different $T^{\prime}_0$. 
To obtain the largest time step
$dt$, we adopt $T^{\prime}_0=3.5$ K. Next, we inspect the resulting mean field quantities with different
$T^{\prime}_0$. \Fig{condens_meant_comp} demonstrates that the mean field quantities converges at $T^{\prime}_0=3.5$ K.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_meanT_compCs}
\end{center}\caption[]{Effect of different $T^{\prime}_0$ on flow field.
We use $\nu=5\times10^{-4}\m^2\s^{-1}$, $\kappa=\nu$ so that
the Prandtl number is $Pr=\nu/\kappa=1$, and $D=\nu$ so that the Schmit number $Sc=\nu/D=1$.
$N_{\rm grid}=16^3$, $N_{\rm s}=N_{\rm grid}/2$. A monodisperse initial distribution of droplets
with $r_{\rm ini}=10\,\mu$m is adopted. The buoyancy force is not invoked.
$\urms$ and $\bar{\epsilon}$ are almost the same for different $T^{\prime}_0$.
When $T^{\prime}_0 \le 3.5$ K, the code crashes. To obtain the largest time step
$dt$, we adopt $T^{\prime}_0=3.5$ K. Cyan line shows the case with $C_d=0$ and
$T^{\prime}_0=0.25$ K. It demonstrates that temperature does not affect the
flow field when buoyancy force is not included.
}
\label{condens_meanT_compCs}
\end{figure*}

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_meanT_comp}
\end{center}\caption[]{Effect of different $T^{\prime}_0$ on mean field quantities.
Same simulations as in \Fig{condens_meanT_compCs}
The mean field quantities converges at $T^{\prime}_0=3.5$ K.
}
\label{condens_meanT_comp}
\end{figure*}

\subsection{Initial configuration}

Saturated or unsaturated initial configuration is important for the evolution
of the size distribution, especially, when collision is involved in.
\Fig{condens_initialCondition_comp} shows the comparison between saturated
and unsaturated initial configuration.
We also compare the size distribution as shown in \Fig{f_condens_initialCondition_comp}. 

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_initialCondition_comp}
\end{center}\caption[]{Different initial configurations.
}
\label{condens_initialCondition_comp}
\end{figure*}

\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{f_condens_initialCondition_comp}
\end{center}\caption[]{Comparison of the size distribution for
different initial configurations. Same simulations as in
\Fig{condens_initialCondition_comp}.
}
\label{f_condens_initialCondition_comp}
\end{figure}

\subsection{Effect of the buoyancy force}

\Fig{condens_buoyancy_comp2} shows the comparison between the cases
with or without the buoyancy force. The buoyancy enhance the
evaporation. This makes the condensaitonal growth slower compared
with the case without buoyancy force as shown in \Fig{f_grav0Re50_cond_buoyancy}.
It is important to have both buoyancy force and gravity included
when studying condnensational growth. Otherwise, the broadening of
the size distribution will be exaggerated as shown in
\Fig{f_grav0Re50_cond_buoyancy}. This could be due to the
interactions among different time scales, like the sedimentation time
scale, turbulence time scale, and the phase transition time scale.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{condens_buoyancy_comp2}
\end{center}\caption[]{With or without the buoyancy force.
Same simulations as in \Fig{f_grav0Re50_cond_coa_comp}.
}
\label{condens_buoyancy_comp2}
\end{figure*}

\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{f_grav0Re50_cond_buoyancy}
\end{center}\caption[]{Size distribution of the same simulations
as in \Fig{condens_buoyancy_comp2}.
}
\label{f_grav0Re50_cond_buoyancy}
\end{figure}

\subsection{Test of $\nabla q_v$}

We also test the initial $\nabla q_v$ and $\nabla T$ as shown in \Fig{gradascalar}.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{gradascalar}
%\includegraphics[width=\textwidth]{gradTT}
\end{center}\caption[]{Different $\nabla q_v$.
Same simulation as in \Fig{condens_meanT_compCs}.
$T^{\prime}_0 \le 3.5$ K
The case with $\nabla q_v=0.1$ crashes after 14000 (0.9 s)
time steps. The value of the gradient should be chosen as the same order of $q_{v0}$. 
}
\label{gradascalar}
\end{figure*}

\section{Verification of the new implementation}

\Fig{moments_conden_comp_newImplementation} demonstrates that
the newly implemented scheme is the same as the original one.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{moments_conden_comp_newImplementation}
\end{center}\caption[]{Verification of the new implementation.
}
\label{moments_conden_comp_newImplementation}
\end{figure*}

\section{Statistical convergence of $N_{\rm s}/N_{\rm grid}$, a revisit}

\Fig{moment_tur_npar_more} shows that $N_{\rm s}/N_{\rm grid}$ converges
even at 0.05.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{moment_tur_npar_more}
\end{center}\caption[]{3-D turbulence. Convergence study of
$N_{\rm p}/N_{\rm grid}$ based on Run 2C.
Here $N_{\rm grid}=128^3$ is fixed. $n_0=10^{10}\m^{-3}$.
$L=0.25$ m. $Re_{\lambda}=100$. Lognormal distribution
with $r_{\rm ini}=10\um$ and $\sigma=0.2$.
$N_{\rm p}/N_{\rm grid}$ even converges at $0.05$.
}
\label{moment_tur_npar_more}
\end{figure*}
%
\section{$\bar{\epsilon}$ dependency on the spatial resolution}

\Tab{tab:resolution} shows that $\bar{\epsilon}$ depends on
the buoyancy force $B$ and spatial resolution.

\begin{table*}[t!]
\caption{Summary of the simulations.}
\centering
\setlength{\tabcolsep}{3pt}
\begin{tabular}{lccccccccccr}
	Run &  $f_0$ & L (\rm{m}) & $N_{\rm grid}$ & $\nu\,(\m^2\s^{-1})$& $B$& $u_{\rm rms}$ ($\rm{m} \, \rm{s}^{-1}$) & $\mbox{\rm Re}_{\lambda}$ & $\bar{\epsilon}$ ($\rm{m}^2\rm{s}^{-3}$)& $\eta\cdot 10^{-4}$ ($\rm{m}$) & $\tau_{\eta}$ ($\rm{s}$) & $\tau_{\rm eddy}$ ($\rm{s}$) \\
\hline
	A & $0.02$  & 0.5  & $512^3$ & $1.5\times10^{-5}$ & Yes &0.32 &142  & 0.059  & 4.9  & 0.016 & 0.50 \\ %3D/SW512condensP_coag0_grav_buoyancy_n2p5e8_L0p5
B & $0.02$  & 0.5  & $512^3$ & $1.0\times10^{-5}$      & no  &0.27 &158 & 0.036 & 4  & 0.017 & 0.6 \\ %3D/SW512condens0_coag_grav0_turb_shima_f2em2_ngp_a24; ; SW1024condens0_coag_grav0_turb_shima_f2em2_ngp_a24
C & $0.02$  & 0.5  & $80^3$ & $1.5\times10^{-5}$& no &0.27 &196 & 0.015  & 6.9  & 0.032 & 0.60 \\ %3D/SW80condensP_coag0_grav_buoyancy_n2p5e8_f0p02_L0p5_noParticle
\hline
\multicolumn{12}{p{0.85\textwidth}}{
Here, $f_0$ is the amplitude of the random forcing (see text) and $L$ is the domain size.
}
\end{tabular}
\label{tab:resolution}
\end{table*}

\section{Re-check of supersaturation fluctuations}

Since the supersaturation fluctuations are extraordinarily large,
we investigate why it is the case. \Fig{debug_G} shows comparison of
the simulation results with different $G$. \Fig{saturationPressure} shows the saturated
vapor pressure $e_s$ and the saturated water vapor $q_{vs}$.

\begin{figure}[t!]\begin{center}
\includegraphics[width=0.5\textwidth]{debug_G}
\end{center}\caption[]{comparison of
the simulation results with different $G$.
$G=1.17\times10^{-10}\m^2\s^{-1}$ is the value used in \cite{Gotzfried17} and
$G=1.2\times10^{-12}\m^2\s^{-1}$ is the value used in \cite{saito2017turbulence}.
}
\label{debug_G}
\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=0.5\textwidth]{saturationPressure}
\end{center}\caption[]{
$G=1.2\times10^{-12}\m^2\s^{-1}$.
Saturated vapor pressure $e_s$ and the saturated water vapor $q_{vs}$.
%3D/SW16condensP_coag0_buoyancyMean
}
\label{saturationPressure}
\end{figure}

Temperature fluctuations is determined by the $Pr$ as shown in
\Fig{kappa_comp}. When $Pr=0.1$, $T_{\rm rms}$ and $s_{\rm rms}$ are
close to the observed values \cite{siebert2017supersaturation}.
The large $s_{\rm rms}$ is due to large $T_{\rm rms}$ as can be seen
from \Fig{ttc_comp}, where a simulation with only the passive scalar
equation of $T$ and the N-S equation is conducted.

\begin{figure}[t!]\begin{center}
\includegraphics[width=0.5\textwidth]{kappa_comp}
\end{center}\caption[]{Different $Pr$.
$G=1.2\times10^{-12}\m^2\s^{-1}$.
}
\label{kappa_comp}
\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=0.5\textwidth]{ttc_comp}
\includegraphics[width=0.5\textwidth]{timeSeries_ttc}
\end{center}\caption[]{Only the passive scalar equation of $T$ is solved.
}
\label{ttc_comp}
\end{figure}

\section{Initial $\nabla T$ and $\nabla q_v$}

\Fig{timeSeries_gradQv} shows time series for
different $\nabla q_v (t=0)$.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{timeSeries_gradQv}
\end{center}\caption[]{Different $\nabla q_v(t=0)$ when $\nabla T(t=0)=1\,K\rm{m}^{-1}$.
}
\label{timeSeries_gradQv}
\end{figure}

\Fig{timeSeries_gradT} shows time series for
different $\nabla T (t=0)$.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{timeSeries_gradT}
\end{center}\caption[]{Different $\nabla T(t=0)$ when $\nabla q_v(t=0)=10^{-4}$.
}
\label{timeSeries_gradT}
\end{figure}

\section{Volume average}

Since updraft cooling is not included in the model,
either adopting $\left<T\right>$, $\left<q_v\right>$
or $T_{\rm env}$, $q_{v,\rm env}$ should be identical.
As we shown in \Fig{timeSeries_volumeMean}, it is indeed the case.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{timeSeries_volumeMean}
\end{center}\caption[]{Comparison of the time series for different
ways of calculating $B$.
}
\label{timeSeries_volumeMean}
\end{figure}

We check the calculation of $\left<T\right>$ (volume average)
in \Fig{volumeMean_debug}.

\begin{figure}[t!]\begin{center}
\includegraphics[width=0.5\textwidth]{volumeMean_debug}
\end{center}\caption[]{Comparison of $\left<T\right>$
and $\bar{T}$ to check the calculation of $\left<T\right>$.
%nils/particles_coag/3D/SW16condensP_coag0_buoyancy_debug
}
\label{volumeMean_debug}
\end{figure}

\section{Removed from paper}

The relative standard deviation of $\sigma_A$ is shown in
\Fig{sigma_areaNorm_f_comp}.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{sigma_areaNorm_f_comp_gradT1_gradQv1em4}
\end{center}\caption[]{Time evolution of $\sigma_A$ for different
$\bar{\epsilon}$ and $Re_{\lambda}$.
Same simulations as in \Fig{timeSeries_gravRe50_cond_e_comp}.
} 
\label{sigma_areaNorm_f_comp}
\end{figure}

Time series of the field quantities from simulations with or without buoyancy force
is shown in \Fig{timeSeries_gravRe50_cond_e_comp_buoyancy0}. The buoyancy force
has a minor effect on the condensational growth. This is close to the conclusion
from the result of \cite{Kumar14}.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{timeSeries_gravRe50_cond_e_comp_buoyancy0}
\end{center}\caption[]{Time averaged power spectra of turbulence
for different $\bar{\epsilon}$ and $Re_{\lambda}$.
Same simulations as in \Fig{timeSeries_gravRe50_cond_e_comp}.
} 
\label{timeSeries_gravRe50_cond_e_comp_buoyancy0}
\end{figure}

\begin{table}[t!]
\caption{Summary of the simulations.}
\centering
\setlength{\tabcolsep}{3pt}
\begin{tabular}{lr}
Quantity ($\rm{s}$) & DNS \\
\hline
$\tau_{\eta}$ & 0.16 \\ 
$\tau_{\rm{p}}$ & 0.002 \\ 
$\tau_{\rm{ph}}=(4\pi D^{\prime}n_0r_{\rm ini})^{-1}$ & 9.448 \\ 
$\tau_{\rm{eddy}}$ & 0.246 \\ 
$\tau_{\rm{ev}}=r_{\rm ini}^2/(2Gs_0)$ & 25.433 \\ 
\hline
\multicolumn{2}{p{0.48\textwidth}}{
Here, $n_0=1.25\times10^8\,\rm{m}^{-3}$ and $r_{\rm ini}=10\,\mu$m.
$D^{\prime}=6.74\times10^{-6}\,\rm{m}^2\rm{s}^{-1}$ is the modified vapor diffusivity,
and $G=1.17\times10^{-10}\,\rm{m}^2\rm{s}^{-1}$ (we adopt values from Table 2 of \cite{Gotzfried17}). 
}
%3D/python/initial_time.py
\end{tabular}
\label{Swarm_Rey}
\end{table}

\cite{desai2018influence}
suggested that the variations of the mean radius and the mean
number density contribute significantly to the increase rate
of the size distribution width, not to the mean growth. It is
therefore important to study the source of the variability,
for example, particle clustering due to turbulence and gravity.
The discretized form of the relaxation time can be written as
\begin{equation}
	\tau_s^{-1}=4 \pi \rho_w A_2 A_3 \sum_{k=1}^{r_{\rm max}} f_k r_k\delta r_k. 
\label{eq:phase2}
\end{equation}

To conserve the mass of liquid water, 
we define the mixing ratio of total water 
as $q_{tot}=q_l+q_v+q_a$, where $q_l$ is the mixing 
ratio of liquid water defined as \citep{saito2017turbulence},
\EQ
q_l\left(\xx,t\right)= 
\frac{4\pi\rho_l}{3\rho_a (\Delta x)^3}\sum_{j=1}^{N_{\triangle}}  
r\left(t\right)^3. 
\EN

%XY: removed the following two figures, just discussion
\Fig{DamKohler_gravRe50_cond_e_comp} shows the time
evolution of $Da$.
The amplitude of the corresponding $Da$ calculated from different configurations
varies from 0.05 to 0.15.
$Da$ increases with increasing $Re_{\lambda}$ and decreases with increasing
$\bar{\epsilon}$.
Consistent with the size distribution dependency on  $Re_{\lambda}$,
larger $Da$ results in wider size distribution.
When $Re_{\lambda}$ is fixed, $Da$ decreases with increasing $\bar{\epsilon}$.
The amplitude of $Da$ for the case of $\bar{\epsilon}=0.005\,\rm{m}^2\rm{s}^{-3}$
is larger than the one for the case of $Re_\lambda=130$, which is inconsistent
with the time evolution of droplet size distributions as shown in
\Fig{f_gravRe50_cond_e_comp}.
We attribute this to the fact that entrainment is not included, $Da\ll1$ in our simulations.
Therefore, $Da$ is not representative to characterize the different configurations.
The ratio between $T_L$ and $\tau_{\rm cond}$ better characterize the scale
separation since $\tau_{\rm cond}\sim s_{\rm rms}^{-1}$. 
\Fig{cond_gravRe50_cond_e_comp} shows the time evolution of
the ratio between $T_L$ and $\tau_{\rm cond}$, which 
is consistent with the size distribution dependency on $Re_{\lambda}$ and
$\bar{\epsilon}$.


%GS:Not interesting to see the time evoloution, we expect it to be approx constant. The ratio varies between 0.05 and 0.15 for all experiments, thus tau_phase is much larger than tau_turb 
%XY: removed the following two figures, just discussion
\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{DamKohler_gravRe50_cond_e_comp_gradT1_gradQv1em4}
\end{center}\caption[]{Time evolution of $Da$ for different
$\bar{\epsilon}$ and $Re_{\lambda}$.
Same simulations as in \Fig{ppower_FixE0p03}.
} 
\label{DamKohler_gravRe50_cond_e_comp}
\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{cond_gravRe50_cond_e_comp_gradT1_gradQv1em4}
\end{center}\caption[]{Time evolution of the ratio between $T_L$
and $\tau_{\rm cond}$for different $\bar{\epsilon}$ and $Re_{\lambda}$.
Same simulations as in \Fig{ppower_FixE0p03}.
%GS:again, don't see the point in time series of this.
} 
\label{cond_gravRe50_cond_e_comp}
\end{figure}

\subsection{Diagnostics}
%GS:Not clear why you have what you have here as diagnostics,
%XY: These are diagnostics from the simulation, which are not part of the model.
%XY: I have now added more explanations. 

Condensational growth of cloud droplets involves in phase transition
characterized by the phase transition time scale $\tau_{\rm phase}$ \citep{2011_lamb},
\begin{equation}
\tau_{\rm phase}^{-1}(t) = 4 \pi \rho_l G \int_{0}^{\infty} r f dr,
\label{eq:phase}
\end{equation}
where $f(r,\xx,t)$ is the size distribution of cloud droplets directly from the simulations.
$\tau_{\rm phase}$ is inversely proportional to the mean radius and
mean number density, indicating that larger cloud droplets or droplets
with high number density can ``eat up'' the water vapor quickly so that
the system can restore to a equilibrium state faster.
$\tau_{\rm phase}$ is determined by the first moment
of the size distribution. Similarly, the size distribution
is also affected by $\tau_{\rm phase}$. 
Another time scale related to the
condensational growth is the evaporation/condensation time scale
\citep{siewert2017statistical},
%GS:This is not really another time scale it is not just the integrated version?
%XY: Not exactly, $s_{\rm rms}$ is introduced here.
\begin{equation}
\tau_{\rm cond}(t) = \frac{\bar{r}^2}{2Gs_{\rm rms}},
\label{eq:con}
\end{equation}
which is determined by the mean radius $\bar{r}$ and the
rms supersaturation $s_{\rm rms}$.
The Damk\"{o}hler number is the ratio of the fluid time scale
to the characteristic thermodynamic time scale associated with
the evaporation process
\EQ
%GS: use \tau for time scales
{\rm Da}=\frac{\tau_L}{\tau_{\rm phase}},
\label{Damkohler}
\EN
where $\tau_L$ is the turbulence integral time scale defined as
$\tau_L=2\pi/(k_f u_{\rm rms})$ ($k_f$ is the turbulence forcing scale
and $u_{\rm rms}$ is the turbulence rms velocity).
$Da\ll1$ characterizes the homogeneous mixing and $Da\gg1$ characterizes
the inhomogeneous mixing \citep{beals2015holographic}.
Homogeneous mixing is facilitated by condensation/evaporation
without entrainment of dry air. Inhomogeneous mixing corresponds
to entrainment of dry air, which is more complicated. 
In the present study, we focus on homogeneous mixing.
%GS:Damköhler number is important here but not the way you describe it,
%GS:here it is simply the condesational time scale and the turbulence time sale that matters.
%XY: what about now?

We only considered the homogeneous mixing ($Da \ll 1$) here. The
inhomogeneous mixing will lead to very rapid changes of supersaturation fluctuations,
which may result in fast broadening of the size distribution \citep{Kumar14}. This
will be explored further in the future studies. 

\section{linear check}

\begin{figure}[h!]\begin{center}
\includegraphics[width=0.49\textwidth]{ssatRe}
\end{center}\caption[]{Linear check for $s_{\rm rms}$ ($[0.15\%, 0.26\%, 0.45\%]$)
versus $\rm{Re}_\lambda$ (45, 78, 130).} 
\label{ssatRe}
\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\textwidth]{sigma_Re}
\end{center}\caption[]{Linear check. [12.00, 9.22, 7.50]; [1.49, 3.41, 7.50]} 
\label{sigma_Re}
\end{figure}

\section{\obsb{Numerical time step}}

\obsb{Figure~\ref{timeSeries_gravRe50_cond_dt} shows the time series
from different numerical time steps of Run A in \Tab{Swarm_Rey}.
Different time steps reveal the same results, suggesting that
the time step is determined by the smallest physical time
step in our simulation.}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\textwidth]{timeSeries_gravRe50_cond_dt}
\end{center}\caption[]{Time series of the field quantities:
(a) $T_{\rm rms}$, (b) $q_{v},\rm{rms}$, and (c) $s_{\rm rms}$.
Run A with different time step.
} 
\label{timeSeries_gravRe50_cond_dt}
\end{figure}


\section{\obsb{Acceleration}}

\obsb{Figure~\ref{timeSeries_arms} shows that the buoyancy acceleration
is about two orders of magnitude smaller than the fluid acceleration.}

\begin{figure}[t!]\begin{center}
\includegraphics[width=0.5\textwidth]{timeSeries_arms}
\end{center}\caption[]{Ratio between the fluid acceleration $a_{\rm rms}$
and the buoyancy acceleration $B_{\rm rms}$.
Run A in \Tab{Swarm_Rey}.
%SW128condensP_coag0_grav_buoyancy_n2p5e8_gradT1_gradQv1em4_digit
} 
\label{timeSeries_arms}
\end{figure}

\section{Response}

Also, since gravity is a potential
field, we estimate the fluid acceleration by superpose the gravity acceleration
and the forced turbulence acceleration.


%---------------------------
\begin{acknowledgments}
We acknowledge support from Research Council of Norway FRINATEK
Grant No. 231444.
\end{acknowledgments}

$G\approx10^{-11}\,\rm{m}^2\rm{s}^{-1}$
%\bibliographystyle{unsrt}
%\bibliographystyle{plainnat}
\bibliography{supersat}

\end{document}
