Computer exercise 3: Survival Probabilities of new particles


This is an exercise to calculate "survival probabilities" of newly formed particles based on the approach described by Westervelt et al. (2013).

From Westervelt et al. (2013): "Dynamics of new particle formation from vapor to cloud droplet. Two competing processes determining the fate of freshly formed atmospheric nuclei are condensational growth (eventually forming CCN) or coagulational scavenging, which results in the loss of the nuclei."

New particles that grow beyond a molecular cluster (~3 nm) can continue growing through net condensation of vapor molecules that increase its size. Given sufficient amount of condensing vapor, these particles will grow to a size (~50--100 nm) large enough to seed cloud droplets. However, some of these smaller, growing particles will be scavenged by larger particles (via process of coagulation) that might have already been destined to become cloud dropets. The outcome -- whether these particles are scavenged before fully grown -- can have a significant impact on the cloud droplet number concentration and cloud droplet sizes, which in turn affects precipitation patterns, cloud lifetimes, and the global energy balance.

From Westervelt et al. (2013): Observed nucleation event in Pittsburgh, PA, on 16 April 2002. "Contours represent values of the size distribution function plotted against particle diameter in m (ordinate) and time in UTC (abscissa). Dashed line represents the diameter growth trajectory."

How many 3 nm particles will grow large enough to possibly serve as seeds for cloud droplet formation? ~ 50--100 nm are typical sizes at which this becomes possible; we will use 100 nm. We can characterize this quantity with a number referred to as the "survival probability" ($\Pr$). If, for instance, if $J_3$ denotes the number of new particles (3 nm size) formed per volume per time, the survival probability can tell us the rate of appearance of $x$ nm sized particles:

$$ J_{x} = \Pr_{3 \to x} \, J_{3} $$

Defining individual transition probabilities from size bin $k$ to $k+1$. How can we calculate this probability? We can write the rate of loss of particles [of size $D_{p,k}$ and number concentration $n_k(t) = n(D_{p,k},t)$] to coagulation as a first order process (temporarily using the symbol, $\kappa$ to denote the first-order rate constant, which is the inverse of the characteristic timescale):

\begin{align} L_k^{(\textrm{coag})} = \kappa_k \, n_k(t) = \frac{1}{\tau_k^{(\textrm{coag})}} n_k(t) \end{align}

The fraction of particles of size $D_{p,k}$ remaining after loss by coagulation over time interval $t$ is therefore given by:

\begin{align} \frac{n_k(t)}{n_k(0)} = \exp\left(- \frac{t}{\tau_k^{(\textrm{coag})}} \right) \end{align}

The characteristic time for particles to grow from $D_{p,k}$ to $D_{p,k+1}$ by condensation is denoted as $\tau_{k,k+1}^{(\textrm{cond})}$. The the fraction of particles that remains after loss to coagulation during this period is given by $n_k(t)/n_k(0)$, which can also be recast as a survival probability (or transition probability) $\Pr_{k \to k+1}$:

\begin{align} \Pr_{k \to k+1} = \frac{n_k(\tau_{k,k+1}^{(\textrm{cond})})}{n_k(0)} = \exp\left(- \frac{\tau_{k,k+1}^{(\textrm{cond})}}{\tau_k^{(\textrm{coag})}} \right) \end{align}

Defining timescales of condensation and coagulation. Now let us determine the appropriate expressions for $\tau$, which is generally be given by $n_k / (\partial n_k/\partial t)_{\text{process}}$. The expressions for condensation/evaporation and coagulation are given by:

\begin{align} \left[\frac{\partial n_k(t)}{\partial t}\right]_{\textrm{cond}} &= - \frac{\partial}{\partial D_p}\left(I(D_{p,k}) n_k(t) \right)\\ \left[\frac{\partial n_k(t)}{\partial t}\right]_{\textrm{coag}} &= P_k^{(\textrm{coag})} - L_k^{(\textrm{coag})} = \frac{1}{2} \sum_{j=1}^{k-1} K(D_{p,j},D_{p,k-j})\, n_j\, n_{k-j} - n_k \sum_{j=1}^{\max(k)} K(D_{p,k},D_{p,j})\, n_j \ \ \ \forall k \geq 2 \end{align}

$P$ and $L$ denote production and loss, respectively, and $K$ is the coagulation coefficient. $I$ is the diameter growth rate.

Note that during lecture, $n$ and therefore $I$ were written as a function of volume ($v$, which is $(\pi/6) D_p^3$ for a sphere). In the transition regime, the volumetric growth rate is given by:

$$ I(v) = \frac{dv}{dt} = \frac{2 \pi D_{AB} M_A}{RT \rho_p} D_p\ f(Kn, \alpha) \left(p_{\infty, A} - p_{\textrm{eq}, A}\right) $$

$M_A$ = molecular weight of species A, $D_{AB}$ = diffusion coefficient of $A$ in air (indexed as $B$), $\rho_p$ = particle density, $R$ = gas constant, $T$ = temperature, $p_A$ = vapor pressure of $A$ in equilibrium with the surface ("eq") and in the bulk vapor phase ("$\infty$") far from the surface, and $f(Kn, \alpha)$ is the Fuchs-Sutugin correction factor in air:

$$ \newcommand{\Kn}{\operatorname{Kn}} f(Kn, \alpha) = \frac{0.75\alpha (1+\Kn)}{\Kn^2 + \Kn+0.283\Kn\alpha + 0.75\alpha} \quad \text{where} \quad \Kn = \frac{2\lambda}{D_p} $$

The mean free path of the condensing molecule in air, $\lambda$, can be calculated from $D_{AB}$ but is provided below in the exercise.

$I$ can alternatively be defined as a function of $D_p$ - then referred to as the diameter growth rate:

$$ I(D_p) = \frac{dD_p}{dt} = ? $$

You are asked to derive this expression in Problem 1 below.

Corresponding expressions for characteristic timescales can be calculated as:

\begin{align} \tau_{k,k+1}^{(\textrm{cond})} &= \frac{n_k}{[\partial n_k / \partial t]_{\textrm{cond}}} \approx - \frac{1}{\Delta_{k,k+1} I / \Delta_{k,k+1} D_p}\\ \tau_k^{(\textrm{coag})} &= \frac{n_k}{L_k^{\textrm{coag}}} = \frac{1}{\frac{1}{2} K(D_{p,k},D_{p,k}) n_k + \sum_{j=k+1}^{\max(k)} K(D_{p,k},D_{p,j}) n_j} \end{align} $\Delta_{k,k+1}(x) = x(k+1) - x(k)$ denotes the forward finite difference operator that can be used to numerically approximate a derivative. Note that we only consider losses by coagulation with particles of size $D_{p,k}$ and greater, and the first term in the denominator of the second equation is the rate of loss by self-coagulation (the 1/2 required because of symmetry: $K_{12} = K_{21}$).

Defining the survival probability. Returning back to our original definition, the overall survival probability between any two sizes $D_{p, m}$ to $D_{p, n}$ is the product of all probabilities that span the range:

$$ \Pr_{m \to n} = \prod_{k=m}^{n-1} \Pr_{k \to k+1} = \prod_{k=m}^{n-1} \exp\left( -\frac{\tau_{k,k+1}^{(\textrm{cond})}}{\tau_k^{(\textrm{coag})}} \right) $$

All diameters $\{D_{p,m}, D_{p, m+1}, \ldots, D_{p,n}\}$ to be used for the calculation will be provided in the exercise.

Problem statement

You will estimate aerosol survival probabilities for a set of background concentration of aerosols ($n_k$) and condensing vapor concentrations.

You are provided with a function that returns $D_p$ ($\mu$m), $\Delta D_p$ ($\mu$m), and $n_k$ (cm$^{-3}$) for use with problems 2 and 3. $n_k$ is normalized such that $\sum_k n_k = 1$. This distribution can be multiplied by a value, $N_0$ cm$^{-3}$, such that the total number of particles equals $N_0$ cm$^{-3}$. A function which calculates the collision frequency function (coagulation kernel) for a pair of particle diameters is also provided so that coagulation rates can be easily computed.

Turn in your own 1) report and 2) Octave/MATLAB code which address the following items:

  1. Find the analytical expression for $I(D_p)$ (spherical particles).
  2. Estimate variables $\{\tau_k^{(\textrm{coag})}, \tau_{k,k+1}^{(\textrm{cond})}, \Pr_{k \to k+1}\}$ as a function of environmental parameters (background concentration of particles and vapor pressure driving the condensation). Plot variables as a function of $D_p$ for each scenario described below:
    • a. Fix $N_0$ to 200 cm$^{-3}$. Calculate variables for $\nabla p_A = \{10^{-12}, 10^{-11}, 10^{-10}, 10^{-9}\}$ atm (which is equivalent to $\{1, 10, 100, 1000\}$ ppt). Note that we have adopted the notation $\nabla p_A = p_{\infty, A} - p_{\textrm{eq}, A}$.
    • b. Fix $\nabla p_A$ to $10^{-9}$ atm (1 ppb). Calculate variables for $N_0 = \{50, 200, 500\}$ cm$^{-3}$.
  3. Estimate the overall survival probabilities (a single number for each case) of 3 nm growing to 100 nm. Using results of individual probabilities obtained from Problem 2, explain (in a few sentences) what type of environmental conditions allow new particles to grow to cloud nuclei sizes by discussing their influences on aerosol process.

Environmental conditions and parameters:

  • $T$ = 298 K, $P$ = 1 atm.
  • Properties of condensing vapor (equivalent to that of H$_2$SO$_4)$: diffusion coefficient in air $D_{AB} = 0.1$ cm$^2$ s$^{-1}$, molar mass $M_A = 98.079$ g mol$^{-1}$, mean free path of in air for use in Knudsen number calculation (for transition regime correction) $\lambda = 0.118$ $\mu$m.
  • Particle density is $\rho_p = 1.0$ g cm$^{-3}$.
  • Assume pressure gradient and number size distribution are constant in each scenario for calculation of characteristic timescales.
  • Assume mass accommodation coefficient ($\alpha$) of unity.


  1. Given the range of particle diameters to be examined in this problem, we are in the transition regime for expressing the net rate of migration of molecules to a particle ("molar flow", $J$). $f(Kn, \alpha)$ corresponds to the correction factor to the continuum regime molar flow ($J_c$) necessary for transition regime calculations.
  2. Write an Octave/MATLAB function for $I(D_p)$ in which you can vary the diameter $D_p$ and concentration gradient. The following provides a suggested template (you fill in the {...}):

     %% -----------------------------------------------------------------------------
     function rate = growth_rate(D_p, gradP)
     %% define constants
     %% calculate Fuchs-Sutugin correction factor
     f = {...}
     %% calculate I(Dp)
     rate = {... use f, Dp, etc. ...}
     %% -----------------------------------------------------------------------------

    Note to take care with unit conversions. You can test this function to solution given below. You do not actually have to integrate this rate function to calculate timescales. For the timescale calculation, you are provided $\Delta D_p$ so that you can calculate $D_{p,k+1} = D_{p,k} + \Delta D_{p,k}$. Further hint: $\Delta_{k,k+1} (I) = I(D_{p,k+1}) - I(D_{p,k})$ and $\Delta_{k,k+1} (D_p) = D_{p,k+1} - D_{p,k} = \Delta D_{p,k}$.

  3. Write an Octave/MATLAB function for the coagulation loss as a function of the number size distribution. Suggested template (you fill in the {...}):

     %% -----------------------------------------------------------------------------
     function kappa = coag_loss_coef(n, Dp)
     kmax = length(Dp);
     kappa = zeros(size(Dp));
     for k = 1:kmax
         j = (k+1):kmax
         kappa(k) = {... you will need to call coag_kernel() ...}
     %% -----------------------------------------------------------------------------

    Note to take care with unit conversions. To provide clarification, n used in the function above is n(k)=N_0*ndist(k) using N_0 and ndist as defined in the example code (provided below).

  4. Use these two functions to calculate your $\tau_k$s.
  5. Calculate your $D_p$-dependent probabilities $\Pr_k$s from $\tau_k$s. Do not take the final product of probabilities (let their dependence on $D_p$ remain) until you get to Problem 3.
  6. For interpretation, the main variables to focus on in this exercise is $\nabla p_A$, which drives the rate of condensation, and the background concentration of pre-existing particles $n_k$ (low concentration regimes might be considered "clean" while high concentration regimes might be considered "polluted"), which drives the rate of coagulation. The pressure gradient $\nabla p_A$ for condensation typically results from the rapid reactions producing condensable vapor, which results in a bulk-phase vapor pressure $p_{\infty,A}$ that is greater than the equilibrium vapor pressure over the surface of a particle $p_{\textrm{eq},A}$.
  7. Remember to use element-wise operators in Octave/MATLAB where necessary. A common error is to use / instead of ./.

Functions provided:

  • generate_sizehist.m: returns $D_p$ [$\mu$m], $\Delta D_p$ [$\mu$m], and $n_k$ [cm$^{-3}$]. $n_k$ is normalized such that $\sum_k n_k = 1$. Download here.
  • coag_kernel.m: $K(D_{p,1}, D_{p,2})$ receives diameters in units of [m] and returns coagulation coefficients in units of [m$^3$ s$^{-1}$]. Download here.

Example usage:

[Dp, dDp, ndist] = generate_sizehist;

N_0 = [50, 200, 500]; % [cm^-3]

subplot(1, 2, 1)
plot(Dp, ndist);
set(gca, 'xscale', 'log')
xlabel('D_p [\mum]')
ylabel('Normalized units')
subplot(1, 2, 2)
plot(Dp, ndist * N_0)
set(gca, 'xscale', 'log')
xlabel('D_p [\mum]')
ylabel('Number concentration [cm^{-3}]')

Your growth rate function can be checked against the following solution. Here I have named my function growth_rate, which receives arguments of $D_p$ [m] and $\nabla p_A$ [atm]. All other parameters are defined within the function.

%% integrate over {0, 1, 2, ..., 18000} seconds with an initial diameter of 3 nm 
%%   for pressure gradients of 1e-12 and 1e-10 atmospheres.

%% define values
t = linspace(0, 5) .* 3600; % [s]
Dp0 = 3e-9; % [m]
gradP1 = 1e-12; % [atm] equivalent to 1 ppt
gradP2 = 1e-10; % [atm] equivalent to 100 ppt

if(exist('OCTAVE_VERSION', 'builtin') ~= 0)
  %% integrate (Octave)
  y1 = lsode(@(Dp, t) growth_rate(Dp, gradP1), Dp0, t);
  y2 = lsode(@(Dp, t) growth_rate(Dp, gradP2), Dp0, t);
  %% integrate (MATLAB)
  [~,y1] = ode45(@(t, Dp) growth_rate(Dp, gradP1), t, Dp0);
  [~,y2] = ode45(@(t, Dp) growth_rate(Dp, gradP2), t, Dp0);

%% plot
plot(t / 3600, [y1,y2] * 1e6)
set(gca, 'yscale', 'log')
xlabel('Time [h]')
ylabel('D_p [\mum]')
legend({'\nabla p_A = 10^{-12} atm', '\nabla p_A = 10^{-10} atm'}, 'location', 'northwest')

Check your solutions

For Problem 2a (note x-scale is limited between 1 nm and 100 nm):



Westervelt, D. M., J. R. Pierce, I. Riipinen, W. Trivitayanurak, A. Hamed, M. Kulmala, A. Laaksonen, S. Decesari, and P. J. Adams (2013): “Formation and Growth of Nucleated Particles into Cloud Condensation Nuclei: Model-Measurement Comparison.” Atmospheric Chemistry and Physics, 13 (15): 7645–63. doi:10.5194/acp-13-7645-2013.