Aerosol condensation process and its verification test: impact of different Gauss-Hermite quadrature points on the mass transfer coefficient

$\color{red}{\text{Last update: 07-25-2019, test strategy, result and presentation need revision.}}$

1 Introduction

Condensation represents the process that the vapor of a compound condenses onto the surface of existing particle population. MAM currently considers two trace gases: sulfur acid gas ($\rm H_2SO_4$, and its condensed particles $\rm SO_4$) and a single lumped gas species of secondary organic aerosol (SOAG, and its condensed particles SOA \cite{liu-2012-gmd}).

The key parameter of condensation process is the mass transfer coefficient of species $l$ ($l \in \{ {\rm SOAG, H_2SO_4} \}$) to/from particle population (unit: $\rm s^{-1}$). By assuming that the particles are spherical and the size distribution of the particle population is described by a few log-normal modes, MAM calculates the mass transfer coefficient of species $l$ in mode $i$ ($C_{_{l,i}}$) as \cite{fuchs-1971}:

\begin{equation} \label{kmt} C_{_{l,i}} = \int_{-\infty}^{\infty} 2 \pi D_{_p} \mathbb{D}_{_{g,l}} f(K_{_{n,l}},\alpha_{_{l}}) n_{_{i}}(\ln D_{_{p}}) d\ln D_{_{p}} \,, \end{equation}

where

\begin{itemize} \item $D_{_{p}}$ is the aerosol particle diameter (unit: m); \item $n_{_{i}}(\ln D_{_{p}})$ is the log-normal size distribution of mode $i$; \item $\mathbb{D}_{_{g,l}}$ is the gas diffusivity (unit: $\rm m^2 \cdot s^{-1}$). It is calculated by \cite{fuller-1966-iec}: \begin{equation} \label{gasd} \mathbb{D}_{_{g,l}} = \frac{10^{-7}T^{1.75} \sqrt{\frac{1}{M_{_{w,l}}}+\frac{1}{M_{_{w,air}}}}} {p \left( \mathbb{V}_{_{D,l}}^\frac{1}{3}+ \mathbb{V}_{_{D,air}}^\frac{1}{3} \right)^2} \,, \end{equation} Here $T$ is air temperature (unit: K); $p$ is air pressure (unit: $\rm atm$ (the standard temperature and pressure)); $\mathbb{V}_{_{D,l}}$ (unitless) is the molecular diffusion volume; $M_{_{w,l}}$ is the molecular weight (unit: $\rm g \cdot mol^{-1}$). Both $\mathbb{V}_{_{D,l}}$ and $M_{_{w,l}}$ are constants for a given species $l$; \item $K_{_{n,l}}$ is the Knudsen number (unitless). It is calculated by: \begin{equation} \label{knu} K_{_{n,l}} = \frac{2\lambda_{_{l}}}{D_{_{p}}} \,, \end{equation} where $\lambda_{_{l}}$ is the mean free path of vapor molecules in air (unit: m) and it is calculated by: \begin{equation} \label{mfp} \lambda_{_{l}} = \frac{3\mathbb{D}_{_{g,l}}}{v_{_{l}}} \,, \end{equation} where $v_{_{l}} = \sqrt{\frac{8RT}{\pi M_{_{w,l}}}}$ is the mean speed of vapor molecules (unit: $\rm m \cdot s^{-1}$); \item $\alpha_{_{l}}$ is the accommodation coefficient (unitless), which is a constant for a given species $l$; \item $f(K_{_{n,l}}, \alpha_{_{l}})$ is a correction factor for the transition regime where the mean free path of vapor molecules is close to the particle diameter. It takes the formula below: \begin{equation} \label{fka} f(K_{_{n,l}}, \alpha_{_{l}}) = \frac{0.75 \alpha_{_{l}} (1+K_{_{n,l}})}{0.75 \alpha_{_{l}} + K_{_{n,l}} (1+K_{_{n,l}}+0.283 \alpha_{_{l}})} \,. \end{equation} \end{itemize}

Currently MAM calculates the $C_{_{H_2SO_4,i}}$ through Eq. \eqref{kmt}. Due to the large uncertainty of speciation and physical properties of SOAG/SOA, MAM assumes $C_{_{SOAG,i}} \approx 0.81 C_{_{H_2SO_4,i}}$ for simplicity. More details about the model parameterization, assumption and numerical implementation can be found in a separate model documentation. MAM uses the Gauss-Hermite quadrature rule to evaluate the Eq. \eqref{kmt}. In order to do so, changing variables is involved (the detailed derivations can be found in the appendix of a separate model documentation) and the final formula of Eq. \eqref{kmt} before using the Gauss-Hermite quadrature rule reads:

\begin{align} J(\ln D_{_p}) &= D_{_p} * f(K_{_{n,l}},\alpha_{_{l}}) \,, \\ x &= \frac{\ln D_{_p} - \left( \ln D_{_{gn,w,i}} + \beta \ln \sigma_{_{g,i}}^2 \right)} {\sqrt{2}\ln \sigma_{_{g,i}}} \,, \label{dp2x} \\ C_{_{H_2SO_4,i}} &= 2\sqrt{\pi} \mathbb{D}_{_{g,H_2SO_4}} N_{_{t,i}} \exp \left( \beta \ln D_{_{gn,w,i}} + \frac{\beta^2}{2} \ln \sigma_{_{g,i}}^2 \right) \int_{-\infty}^{\infty} \exp^{-x^2} \frac{J(\sqrt2 x \ln \sigma_{_{g,i}} + \ln D_{_{gn,w,i}} + \beta \ln \sigma_{_{g,i}}^2)}{\exp \left[ \beta \left( \sqrt2 x \sigma_{_{g,i}} + \ln D_{_{gn,w,i}} + \beta \ln \sigma_{_{g,i}}^2 \right) \right]} dx \,, \label{ghq} \\ \end{align}

where $\beta$ is an internally calculated float number between 1 and 2; $D_{_{gn,w,i}}$ and $\sigma_{_{g,i}}$ are the geometric mean wet diameter and geometric standard deviation for mode $i$.

2 Goals

MAM uses 2 Gauss-Hermite quadrature points to evaluate the integral in Eq. \eqref{kmt} but it remains unknown whether 2 quadrature points can provide an accurate enough solution. Thus in this test, we aim at examining the impact of different Gauss-Hermite quadrature points on the $C_{_{l,i}}$ in MAM. Since MAM assumes $C_{_{SOAG,i}} \approx 0.81 C_{_{H_2SO_4,i}}$, we then focus on the verification of $C_{_{H_2SO_4,i}}$.

3 Steps to establish the verification test

This verification test uses a box model that inherits the data structure and condensation subroutines from MAM in E3SM. The box model serves as a driver program to specify the necessary parameters in order to calculate the test and reference solutions.

3.1 Test solution

We use the condensation subroutines in MAM but 2, 4, 10, 15, 20, 30 and 40 Gauss-Hermite quadrature points to evaluate the Eq. \eqref{kmt} separately. These solutions are denoted as $C_{_{TEST}}$. The accuracy of using Gauss-Hermite quadrature rule to evaluate the integral $\int_{-\infty}^{\infty} f(x)e^{-x^2} dx$, where f(x) is a polynomial, has been verified in a separate verification test.

3.2 Reference solution

We use condensation subroutines in MAM and 50 Gauss-Hermite quadrature points to evaluate the Eq. \eqref{kmt}. This solution is denoted as $C_{_{REF}}$.

3.3 Parameters required for the verification test

In this verification test, the following quantities are specified:

\begin{enumerate} \item $\alpha_{_{H_2SO_4}} = 0.65$ \cite{poschl-1998-jpca}. \item Temperature $T$ (unit: K) and pressure $P$ (unit: Pa), for which 5 pairs of values are used to represent the typical conditions from troposphere to stratosphere: \begin{align} \rm T&: 310, 298, 273, 230, 200 \\ \rm P&: 10^5, 8.5 \times 10^4, 7 \times 10^4, 2 \times 10^4, 10^4 \end{align} \item Molecular diffusion volume $\mathbb{V}_{_{D,H_2SO_4}} = 42.88$ and $\mathbb{V}_{_{D,air}} = 20.1$ \item Molecular weight $M_{_{w,H_2SO_4}} = 98.0783997$ and $M_{_{w,air}} = 28.966$. \item Number concentration of particle population within each mode ($N_{_{t,i}}$, unit: $\# \cdot m^{-3}$). We set the values based on our previous experience and they are summarized in the table below. \item Geometric mean dry diameter ($D_{_{gn,d,i}}$, unit: $\mu m$) and geometric standard deviation ($\sigma_{_{g,i}}$) for each mode. The values of $D_{_{gn,d,i}}$ and $\sigma_{_{g,i}}$ are chosen to represent the midrange values for tropospheric aerosol observations, and are listed in the table below \cite{dick-2004-jgr}. The geometric mean wet diameter ($D_{_{gn,w,i}}$) is assumed to be 1.2 times the $D_{_{gn,d,i}}$ in this verification test for simplicity. $D_{gn,w,i}$ is the actual value used to evaluate the integral in Eq. \eqref{kmt}. \item Gauss-Hermite quadrature points ($\mathbb{R_{_i}}$) and their associated weights ($\mathbb{W_{_i}}$). They are generated using the Python function **numpy.polynomial.hermite.hermgauss** with 16 decimal digits \cite{Python-ghq}. \end{enumerate}
$$ \text{Table 1. List of typical geometric mean dry diameter ($D_{_{gn,d,i}}$), geometric standard deviation ($\sigma_{_{g,i}}$) and number concentration of } \\ \text{particle population ($N_{_{t,i}}$) for each log-normal mode.}\\ \begin{array}{c|c|c|c} \hline \text{Mode} & D_{_{gn,d,i}}, \, (\mu m) & \sigma_{_{g,i}} & N_{_{t,i}}, \, (\# \cdot m^{-3}) \\ \hline \text{Aitken} & 0.026 & 1.6 & 10^9 \\ \text{Accumulation} & 0.11 & 1.8 & 10^8 \\ \text{Coarse} & 2.0 & 1.8 & 10^5 \\ \text{Primary Carbon} & 0.05 & 1.6 & 2 \times 10^8\\ \hline \end{array} $$

With these parameters known, totally 20 (= 5(groups of T/P) × 4(modes for aerosols)) verification tests are performed.

3.4 PASS/FAIL criterion

In each test, the driver program calls the condensation subroutines in MAM to calculate the $C_{_{TEST}}$ and $C_{_{REF}}$. Then the relative error between $C_{_{TEST}}$ and $C_{_{REF}}$ can be calculated by:

\begin{equation} \phi = \frac{\left| C_{_{TEST}} - C_{_{REF}} \right|}{C_{_{REF}}} \end{equation}

We further perform a linear regression between number of quadrature points and $\phi$. If the slope is negative for all the 20 tests, meaning the $\phi$ decreases with more quadrature points, an overall PASS is then issued.

4 Source codes and configuration files of verification test

The verification test requires a test driver to perform the verification test and the configuration files for compilation and execution of cmdv-test-runner (e.g., the yaml file). If we go the root directory of CMDV-Verification folder that is cloned from GitHub Repo, the test driver code (driver.F90) is under $\rm \color{blue}{tests/mam/restructed\_tests/F90\_SRC/test\_drivers/condensation/gauss\_hermite\_convergence}$ and the yaml file is under $\rm \color{blue}{tests/mam/restructed\_tests/yaml\_files\_for\_test\_runner/condensation/gauss\_hermite\_convergence.test.yaml}$.

5 Verification results

If the verification test runs successfully, we should obtain the netcdf file named "condensation_verification_test.nc" with the mass transfer coefficient calculated by 2, 4, 10, 15, 20, 30, 40 and 50 Gauss-Hermite quadrature points. Assuming that the netcdf file is put under the "postprocess_input" folder, we could show the plots of $\phi$ against the number of quadrature points.

In [1]:
# load library and modules
import numpy as np
import matplotlib.pyplot as plt
import postprocess_scripts.plot_conden_coeff as pc
#turn interactive plotting off
plt.ioff()

#define the file path and name
nc_file = "./postprocess_input/condensation_verification_test.nc"

#define output file path
outpath = "./postprocess_output/"
In [2]:
#make the plot and save the figure (done in the python script)
pc.plot_coeff(nc_file, outpath, ylim=[1.e-18,1], left=0.1, right=1.2, 
              top=0.45, bottom=0.1, hspace=0.4, wspace=0.3, nrows=1,
              ncols=5, fontheight=4.5, loc='upper right')

If a user uses the same settings in Section 3.3 and does not change the source codes, the whole output would look like:

In [3]:
pc.show_figures("./trusted_output/",left=0.1,right=2.5,top=0.9,bottom=0.1,hspace=0.1,wspace=0.1)

Based on the trusted output above, $\phi$ decreases when using more quadrature points. It also seems that in some situations, using fewer quadrature points (the specfic number may differ among modes) could achieve very close result to the reference solution (the curves go straight down). Further increase of quadrature points would not help reduce $\phi$, though $\phi$ is generally smaller than $10^{-15}$. We also observe that using 2 Gauss-Hermite quadrature points leads to $\phi$ between $10^{-4}$ and $10^{-2}$, which is much higher than that of using 10 or more Gauss-Hermite quadrature points. To understand this, we first denote $G(x) = \frac{J(\sqrt2 x \ln \sigma_{_{g,i}} + \ln D_{_{gn,w,i}} + \beta \ln \sigma_{_{g,i}}^2)}{\exp \left[ \beta \left( \sqrt2 x \sigma_{_{g,i}} + \ln D_{_{gn,w,i}} + \beta \ln \sigma_{_{g,i}}^2 \right) \right]}$ in the Eq. \eqref{ghq}. Then we choose a series of $D_{_p}$ and plot the curves of $G(x)$ against $x$ (x is converted from $D_{_p}$ based on the Eq. \eqref{dp2x}). In addition, we fit four polynomials with user-specified degrees (say $\it{porder}$) to the pairs of points $(x, G(x))$. The results look like:

In [4]:
# this cell does the polynomial fit to the G(x)
pc.fit_poly(nc_file, outpath, case_idx=0, porder=[1,3,7,12], 
            xvalues=np.log(np.float_power(10,np.linspace(-9,-4,21))), 
            colors=['forestgreen','steelblue','blueviolet','orange'], 
            xlim=[-14,14], ylim=[1.e-2, 1.e7], 
            xtick=np.arange(-14,21,7), ytick=np.float_power(10,np.linspace(-2,7,10)), 
            yticklabel=[r'$10^{-2}$',r'$10^{-1}$',r'$10^0$',r'$10^1$',r'$10^2$',r'$10^3$', \
                        r'$10^4$',r'$10^5$',r'$10^6$',r'$10^7$'], 
            left=0.1, right=0.9, top=0.35, bottom=0.1, hspace=0.7, wspace=0.3, loc=[0.05,0.05])
In [5]:
# choose a different case_idx for different T/P
pc.fit_poly(nc_file, outpath, case_idx=4, porder=[1,3,7,12], 
            xvalues=np.log(np.float_power(10,np.linspace(-9,-4,21))), 
            colors=['forestgreen','steelblue','blueviolet','orange'], 
            xlim=[-14,14], ylim=[1.e1, 1.e6], 
            xtick=np.arange(-14,21,7), ytick=np.float_power(10,np.linspace(1,6,6)), 
            yticklabel=[r'$10^1$',r'$10^2$',r'$10^3$',r'$10^4$',r'$10^5$',r'$10^6$'], 
            left=0.1, right=0.9, top=0.35, bottom=0.1, hspace=0.7, wspace=0.3, loc=[0.05,0.05])

The results above show that the polynomial does not fit the pairs of points $(x, G(x))$ well until the degree is 12 or higher. This conclusion holds for different temperature and pressure. Recall that using $K$ Gauss-Hermite quadrature points could solve the Eq. \eqref{ghq} exactly when $G(x)$ is a polynomial up to degree $2K-1$. Therefore, it is understandable that using 2 or 4 quadrature points still leads to relatively larger $\phi$ while using 10 or more quadrature points could give a smaller $\phi$ (< $10^{-8}$).

Finally we check whether a PASS/FAIL should be issued based on the criterion in Section 3.4.

In [6]:
# run this cell to issue the PASS/FAIL
pc.check_rate(nc_file)

PASS

References

[1] Liu X., Easter R. C., Ghan S. J. et al., ``Toward a minimal representation of aerosols in climate models: description and evaluation in the Community Atmosphere Model CAM5'', Geoscientific Model Development, vol. 5, number 3, pp. 709--739, 2012. online

[2] NA Fuchs and AG Sutugin, ``High dispersed aerosols: in: Topics in current aerosol research, edited by: Hidy, GM and Brock, JR, vol. 2'', 1971.

[3] Fuller Edward N, Schettler Paul D and Giddings J Calvin, ``New method for prediction of binary gas-phase diffusion coefficients'', Industrial \& Engineering Chemistry, vol. 58, number 5, pp. 18--27, 1966.

[4] P{\"o}schl Ulrich, Canagaratna Manjula, Jayne John T et al., ``Mass accommodation coefficient of H2SO4 vapor on aqueous sulfuric acid surfaces and gaseous diffusion coefficient of H2SO4 in N2/H2O'', The Journal of Physical Chemistry A, vol. 102, number 49, pp. 10082--10089, 1998.

[5] Easter Richard C., Ghan Steven J., Zhang Yang et al., ``MIRAGE: Model description and evaluation of aerosols and trace gases'', Journal of Geophysical Research: Atmospheres, vol. 109, number D20, pp. , 2004.

[6] , ``Python function to generate the Gauss–Hermite quadrature points'', . online