MAM Coagulation Verification Test¶
Component Being Tested¶
Atmosphere, Modal Aerosol Model physics, Coagulation time-stepping
Description¶
Aerosol Modes¶
The Modal Aerosol Model (MAM) tracks aerosol particles via modes, so that a single distribution of particles over a range of diameters can be described by just three terms: a total number concentration, a median particle diameter, and a standard deviation. The distribution, $n$ of particles of diameter $D_i$, for the $i$th mode, is given by
$$ n(\ln D_i) = \frac{N_{\text{t},i}}{\sqrt{2\pi}\ln\sigma_{\text{g},i}} \ \text{exp}\left[-\frac{1}{2}\frac{(\ln D_i -\ln D_{\text{gn},i})^2}{\ln^2 \sigma_{\text{g},i}}\right], $$where $N_{\text{t},i}$ is the total number concentration, $\sigma_{\text{g},i}$ is the standard deviation, and $D_{\text{gn},i}$ is the median particle diameter. The subscript $\text{g}$ denotes a geometric quantity, $i$ denotes the $i$th mode, and MAM tracks three modes. Note that this equation is a little different than a standard normal distribution, in that the random variable is the $\ln D_i$ rather than $D_i$ itself. Similarly, the effective mean is $\ln D_{\text{gn},i}$ and the effective standard deviation is $\ln \sigma_{\text{g},i}$.
As an example of a single aerosol mode, for a total number concentration $N_{\text{t},i} = 10^8$, a median diameter $D_{\text{gn},i} = 0.2 \times 10^{-6}$m, and a standard deviation $\sigma_{\text{g},i} = 1.8$m, the mode would look like this:
# Import packages:
import numpy
import matplotlib.pyplot as plt
import mam_util
from ResultReporter import ResultWriter, ResultReader
from imp import reload # This is python 2/3 compatible
%matplotlib inline
mam_util.plot_single_mode(magnitude=1e8, median=0.2e-6, std_dev=1.8, label="Example")
The three modes tracked by MAM are referred to as the Aitken mode, the Accumulation mode and the Primary Carbon mode. Each mode has its own fixed standard deviation, but the total number concentration and median particle diameters can vary over time. A composite of these modes at a typical distribution might look like this:
aitken = (1e9, 0.04e-6, 1.6, "Aitken")
accum = (1e8, 0.20e-6, 1.8, "Accumulation")
pca = (2e8, 0.08e-6, 1.6, "Primary Carbon")
mam_util.plot_three_modes(aitken, accum, pca)
Commonly in MAM, a normalized number distribution is tracked:
$$ n_{\text{norm},i} = \frac{n(\ln D_i)}{N_{\text{t},i}}. $$Coagulation Rate Coefficients¶
Depending on the Knudsen number, different coagulation kernels apply. The Knudsen number is a dimensionless parameter defined by
$$ \mbox{Kn} = \frac{2\lambda}{D}, $$where $\lambda$ is the mean free path of the fluid surrounding a particle, and $D=\sqrt[3]{D_1^3+D_2^3}$ is the effective diameter of the particle created when two particles of diameter $D_1$ and $D_2$ collide.
$$ \begin{eqnarray} \beta_{nc}(D_1,D_2) &=& 2 \pi\left[\frac{k_\text{B} T C_c}{3\pi\mu D_1} + \frac{k_\text{B} T C_c}{3\pi\mu D_2}\right] (D_1 + D_2), \\ \tilde{\beta}_\text{fm}(D_1,D_2) &=& \left[\frac{3k_\text{B} T}{\rho_\text{p}}\right]^{1/2} \left[D_1^{0.5} + 2\frac{D_2}{D_1^{0.5}} + \frac{D_2^2}{D_1^{1.5}} + \frac{D_1^2}{D_2^{1.5}} + 2\frac{D_1}{D_2^{0.5}}+D_2^{0.5}\right] \end{eqnarray} $$where $C_c$ is the slip correction factor, $\mu$ is the dynamic viscosity of air, $T$ is temperature (K), and $k_\text{B}$ is the Boltzmann constant. $\beta_\text{nc}(D_1,D_2)$ applies in the near continuum regime $(\mbox{Kn} \le 1)$, and $\beta_\text{fm}(D_1,D_2)$ applies in the free-molecular regime $(\mbox{Kn} > 10)$. For each of these, coagulation rate coefficients for intramodal and intermodal coagulation are calculated from
$$ \begin{eqnarray} B_{0,i,i} &=& \left[\frac{1}{2}\int_0^\infty n_{\text{norm},i} \mbox{d}(D_1) \int_0^\infty \beta(D_1,D_2) n_{\text{norm},i} \mbox{d}(D_2)\right] \\ B_{0,i,j} &=& \left[\int_0^\infty n_{\text{norm},i} \mbox{d}(D_{1}) \int_0^\infty \beta(D_1,D_2) n_{\text{norm},j} \mbox{d}(D_2)\right] \\ B_{3,i,j} &=& \left[\int_0^\infty D_1^3 n_{\text{norm},i} \mbox{d}(D_{1}) \int_0^\infty \beta(D_1,D_2) n_{\text{norm},j} \mbox{d}(D_2)\right] \left[\frac{1}{D_{\text{gn},i}^3 \exp\left[\frac{9}{2} \ln^2\sigma_{\text{g},i}\right]}\right] \end{eqnarray} $$where the subscripts $0$ and $3$ correspond to zeroth and third moments. Once $B_{k,\text{nc}}$ and $B_{k,\text{fm}}$ are both known, the final coagulation rate coefficient used by MAM is given by the harmonic mean:
$$ B_k = \frac{B_{k,\text{fm}} B_{k,\text{nc}}}{B_{k,\text{fm}} + B_{k,\text{nc}}} $$Fuchs Kernel¶
In the verification test, the MAM coagulation rate coefficient is compared to the coefficient found using the Fuchs kernel:
$$ \beta\left(D_1,D_2\right) = \frac{2\pi(\textbf{D}_1 + \textbf{D}_2)(D_1 + D_2)}{\left(\frac{D_1 + D_2}{D_1 + D_2 + 2(g_1^2 + g_2^2)^{1/2}} + \frac{8(\textbf{D}_1 + \textbf{D}_2)}{(\bar{c}_1^2 + \bar{c}_2^2)^{1/2}(D_1 + D_2)}\right)}, $$where
$$ \begin{eqnarray} \bar{c}_i &=& \left(\frac{8kT}{\pi m_i}\right)^{1/2} \mbox{m/s} \\ \ell_i &=& \frac{8\textbf{D}_{i}}{\pi \bar{c}_{i}}\mbox{m} \\ g_i &=& \frac{1}{3D_i \ell_i}\left[(D_i + \ell_i)^3 - (D_i^2 + \ell_i^2)^{3/2}\right] - D_{pi} \mbox{m} \\ \textbf{D}_i &=&\frac{kTC_{c_i}}{3\pi \mu D_i} \mbox{m}^2/\mbox{s} \\ C_{c_i} &=& 1+\frac{2\lambda}{D_{pi}}\left[1.257 + 0.4\exp\left(-\frac{1.1D_i}{2\lambda}\right)\right] \end{eqnarray} $$Greg: what is $D_{pi}$?
Verification Test¶
Describe the verification test here.
Variable Definitions¶
| Name | Definition |
|---|---|
| $C_{c_i}$ | slip correction factor for mode $i$ |
| $D_i$ | Particle diameter for mode $i$. The PDF random variable is $(\ln D_i)$ |
| $D_{\text{gn},i}$ | Mean particle diameter for mode $i$ |
| $D_{pi}$ | Unknown |
| $g$ | A subscript denoting geometric quantities |
| $i$ | Mode number |
| $k_\text{B}$ | Boltzmann constant |
| $\mu$ | dynamic viscosity of air |
| $N_{\text{t},i}$ | Total number concentration |
| $n$ | The number distribution of particles of diameter $D_i$ |
| $n_{\text{norm},i}$ | The normailzed number distribution of particles of diameter $D_i$ |
| $\rho_\text{p}$ | Unknown |
| $\sigma_{\text{g},i}$ | Standard deviation for mode $i$ |
| $T$ | temperature (K) |
Discussion¶
Many of the technical details involved in calculating the MAM coagulation rate coefficients and the tests utilized here can be found in the appendix of [1].
Initial Conditions¶
The verification test uses ten different sets of initial number concentrations for the three modes. The sets of initial conditions are described in the following table.
| IC | Accumulation | Aitken | Primary Carbon |
|---|---|---|---|
| 1 | 1.0e+08 | 1.0e+09 | 2.0e+08 |
| 2 | 1.0e+11 | 1.0e+09 | 2.0e+08 |
| 3 | 1.0e+08 | 1.0e+12 | 2.0e+08 |
| 4 | 1.0e+08 | 1.0e+09 | 2.0e+11 |
| 5 | 1.0e+11 | 1.0e+12 | 2.0e+11 |
| 6 | 1.0e+11 | 0.0e+00 | 0.0e+00 |
| 7 | 0.0e+00 | 1.0e+09 | 2.0e+08 |
| 8 | 1.0e+08 | 0.0e+00 | 2.0e+08 |
| 9 | 1.0e+08 | 1.0e+09 | 0.0e+00 |
| 10 | 0.0e+00 | 1.0e+09 | 0.0e+00 |
Results¶
This test is set up to build and run the standalone test driver dd.x compiled from source code coag_b1_driver.F90. It creates output file coag_delNum_delMass.out, and this notebook is responsible for postprocessing and analyzing the results.
#Parse output file in preparation for plotting:
reload(mam_util)
result = mam_util.parse_output_file("coag_delNum_delMass.out")
#First plotting example from verification test 2 (MAM vs RK4)
reload(mam_util)
#Set initial condition (1-10):
IC = 5
#Initialize figure and axis:
figure_size = (8,5)
fig, ax = plt.subplots(1, 1, figsize=figure_size)
#Plot a quantity from the verification solution:
mam_util.plot_quantity_test2(ax=ax, result=result, IC=IC, key="qa(VER,SO4)", label="VER")
#Plot the same quantity from the MAM solution:
mam_util.plot_quantity_test2(ax=ax, result=result, IC=IC, key="qa(MAM,SO4)", label="MAM")
plt.show()
# Second plotting example from verification test 2 (MAM time-split vs RK4)
reload(mam_util)
from ResultReporter import ResultWriter
rw = ResultWriter("coag_results.log")
conv_rates = numpy.zeros((10,12))
# Plot the error (comparing dt=225,450,900,1800 to dt=1):
for IC in range(1,11):
conv_rates[IC-1,:] = mam_util.plot_errors_test2(figure_size, result, IC)
for i in range(conv_rates.shape[1]):
test_name = "IC = %d, %s" % (IC, result["ordered errors"][i])
if conv_rates[IC-1,i] == 0.0:
# This happens when errors are machine zero
test_name += ", TINY errors"
rw.report_test_passed(test_name)
else:
# Convergence rate is finite, check that it is sufficiently large
test_name += ", slope = %g" % conv_rates[IC-1,i]
rw.report_test(test_name, conv_rates[IC-1,i] > 0.8)
rw.finished()
References¶
1. Whitby, E.; McMurry, P.; Shankar, U.; Binkowski, F. Modal aerosol dynamics modeling. 1991.