Aerosol condensation process and its verification test: self-convergence rate of semi-implicit solutions of SOAG/SOA

$\color{red}{\text{Last update: 07-29-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}).

By assuming that the particles are spherical and the size distribution of the particle population is described by a few log-normal modes, the governing equations for the condensation of SOAG read:

\begin{align} C_{_{H_2SO_4,i}} &= \int_{-\infty}^{\infty} 2 \pi D_{_p} \mathbb{D}_{_{g,H_2SO_4}} f(K_{_{n,H_2SO_4}},\alpha_{_{H_2SO_4}}) n_{_{i}}(\ln D_{_{p}}) d\ln D_{_{p}} \,, \label{kmt-h2so4} \\ C_{_{SOAG,i}} &\approx 0.81C_{_{H_2SO_4,i}} \,, \label{kmt-soag} \\ \frac{dq_{_{v,SOAG}}}{dt} &= \sum_{i=1}^{I} C_{_{SOAG,i}} (q_{_{SOAG}}-q_{_{v,SOAG,sat,i}}) \,, \label{qgsoag} \\ \frac{dq_{_{m,SOA,i}}}{dt} &= C_{_{SOAG,i}} (q_{_{v,SOAG}}-q_{_{v,SOAG,sat,i}}) \,, i = 1,2,...,I \,, \label{qasoa} \end{align}

where the definitions of variables in Eq. \eqref{kmt-h2so4} can be found in a separate model documentation; $q_{_{v,SOAG}}$ is the mass mixing ratio of SOAG (unit: kmol of SOAG per kmol of dry air); $q_{_{m,SOA,i}}$ is the mass mixing ratio of SOA in mode $i$; $q_{_{v,SOAG,sat,i}}$ is the saturation vapor mixing ratio of SOAG for mode $i$, considering the existence of oxygenated primary organic aerosol (OPOA).

To solve for $q_{_{v,SOAG}}$ and $q_{_{m,SOA,i}}$, MAM assumes that $C_{_{SOAG,i}}$ are fixed within one time step $\Delta t_{_{cond}}$ ($\Delta t_{_{cond}} = t_{_1} - t_{_0}$). A semi-implicit method is then used to solve Eq. \eqref{qgsoag} and \eqref{qasoa}. More details about the basic notations, model parameterization, assumptions and numerical implementation can be found in a separate ''MAM_Basics'' and model documentation.

2 Goals

It is necessary to check the semi-implicit method for SOAG/SOA is implemented correctly and one criterion is that the relative errors of solutions decrease linearly with the decrease of $\Delta t_{_{cond}}$ in a log-log scale plot (e.g., $1^{st}$-order convergence rate). Thus in this test, we aim at examining the convergence rates of semi-implicit solutions of SOAG/SOA in MAM. We would perform a self-convergence test to facilitate this purpose. Note that for this test, $C_{_{SOAG,i}}$ are fixed over the whole simulation time and we do not consider the existence of $\rm H_2SO_4/SO_4$.

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 MAM and reference solutions.

3.1 MAM solution

According to the mass conservation law, the molar mixing ratio of SOAG/SOA is conserved over any time interval during the condensation process. Thus we could define a constant $q_{_{tot}}$ as:

\begin{equation} q_{_{tot}} = q_{_{v,SOAG}}(t_{_0}) + \sum\limits_{i=1}^{I} q_{_{m,SOA,i}}(t_{_0}) \,. \end{equation}

MAM implemented a dynamical sub-stepping within one $\Delta t_{_{cond}}$. Denote the varying-length sub-steps as $\Delta t_{_{SOAG}}$ ($\Delta t_{_{SOAG}} \le \Delta t_{_{cond}}$, $\Delta t_{_{SOAG}} = t_{_{new}} - t_{_{old}}$, $t_{_0} \le t_{_{old}} \lt t_{_{new}} \le t_{_1}$), the semi-implicit solutions of Eq. \eqref{qgsoag} and \eqref{qasoa} from $t_{_{old}}$ to $t_{_{new}}$ read:

\begin{align} q_{_{v,SOAG}}(t_{_{new}}) &= \frac{q_{_{tot}} - \sum\limits_{i=1}^{I} \frac{q_{_{m,SOA,i}}(t_{_{old}})}{\chi_{_i}^*}} {1 + \sum\limits_{i=1}^{I} \frac{\psi_{_i}^*}{\chi_{_i}^*}} \,, \label{sol-soag} \\ q_{_{m,SOA,i}}(t_{_{new}}) &= \frac{q_{_{m,SOA,i}}(t_{_{old}}) + \psi_{_i}^* q_{_{v,SOAG}}(t_{_{new}})}{\chi_{_i}^*} \,, \label{sol-soa} \end{align}

where $\psi_{_i}^* = C_{_{SOAG,i}} (t_{_0}) \Delta t_{_{SOAG}}$, $\chi_{_i}^* = 1 + \psi_{_i}^* \theta_{_i}^*$ and $\theta_{_i}^* = \frac{q_{_{v,SOAG,equ}} (t_{_0})}{\max \left( q_{_{m,SOA,i}}^*+q_{_{m,OPOA,i}} (t_{_0}), \, 10^{-20} \right)}$. $q_{_{v,SOAG,equ}}(t_{_0})$ is the saturation vapor mixing ratio of SOAG without considering the existence of OPOA. After a few sub-steps, $t_{_{new}} = t_{_1}$ and MAM solves for the solutions of SOAG/SOA after one $\Delta t_{_{cond}}$. More details can be found in a separate model documentation.

3.2 Reference solution

The reference solution is generated by the same method mentioned in Section 3.1 but with $\Delta t_{_{cond}} = 0.01$s. The reason of using $\Delta t_{_{cond}} = 0.01$s to generate the reference solution is that current MAM subroutines do not calculate the condensation of SOAG/SOA if $\Delta t_{_{cond}} \le 10^{-3}$s.

3.3 Parameters required for this verification test

To perform this verification test, the following quantities are specified: %

\begin{enumerate} \item All the parameters used to calculate the $C_{_{H_2SO_4,i}}$. These include: \begin{itemize} \item Accommodation coefficient $\alpha_{_{H_2SO_4}} = 0.65$ (unitless, \cite{poschl-1998-jpca}). \item Temperature $T$, for which we use 273 K as an example. \item Pressure $P$, for which we use 1000 hPa as an exmaple. \item Molecular diffusion volume (unitless) $\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 in mode $i$ ($N_{_{t,i}}$, unit: $\# \cdot m^{-3}$). We set the typical numbers in reality. \begin{equation} N_{_{t,1}} = 10^8, N_{_{t,2}} = 10^9, N_{_{t,3}} = 10^5, N_{_{t,4}} = 2 \times 10^8. \nonumber \end{equation} \item Geometric mean dry diameter ($D_{_{gn,d,i}}$) 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 \cite{dick-2004-jgr, liu-2012-gmd}, and are listed in the ''MAM_Basics'' documentation. The geometric mean wet diameter ($D_{_{gn,w,i}}$) is the actual value used to calculate the $C_{_{l,i}}$ and is assumed to be 1.2 times $D_{_{gn,d,i}}$ for simplicity. \item Gauss-Hermite quadrature points and their associated weights. For this test, both MAM and reference solutions use 2 Gauss-Hermite quadrature points. They are generated using the Python function **numpy.polynomial.hermite.hermgauss** with 16 decimal digits \cite{Python-ghq}. \end{itemize} With $C_{_{H_2SO_4,i}}$ known, $C_{_{SOAG,i}}$ is then approximated by Eq. \eqref{kmt-soag}. \item Simulation length $\Delta t_{_{tot}}$, which is set to 1800 seconds. \item Time step $\Delta t_{_{cond}}$ (unit: second). In this test, we use different $\Delta t_{_{cond}}$ to generate the MAM solutions: \begin{equation} \Delta t_{_{cond}} \in [1800.,900.,600.,300.,100.,60.,30.,15.,5.,1.,0.1,0.01] \nonumber \end{equation} We use $\Delta t_{_{cond}} = 0.01$s to generate the reference solution. \item Density (unit: $\rm kg \cdot m^{-3}$) and molecular weight (unit: $\rm g \cdot mol^{-1}$) of aerosol species ($\rho_{_{l}}$ and $M_{_{w,l}}$), for which the following values are used: $$[\rho_{_{SO_4}}:1770 \ ; \rho_{_{POA}}:1000 \ ; \rho_{_{SOA}}:1000 \ ; \rho_{_{BC}}:1700 \ ; \rho_{_{SS}}:1900 \ ; \rho_{_{DST}}:2600 \ ; \rho_{_{MOA}}:1601] \,, \nonumber$$ $$[M_{_{w,SO_4}}:115 \ ; M_{_{w,POA}}:12 \ ; M_{_{w,SOA}}:12 \ ; M_{_{w,BC}}:12 \ ; M_{_{w,SS}}:58.5 \ ; M_{_{w,DST}}:135 \ ; M_{_{w,MOA}}:250092] \,. \nonumber$$ \item Initial mass fractions of aerosol species in mode $i$ ($f_{_{m,l,i}}$), whose values are shown below (numbers in the parentheses) as the example setting. \begin{align} \text{ Accumulation}:& &\text{SS(50%)},& &\text{SOA(50%)}, \nonumber \\ \text{ Aitken}:& &\text{SS(50%)},& &\text{SOA(50%)}, \nonumber \\ \text{ Coarse}:& &\text{DST(50%)},& &\text{SOA(50%)}, \nonumber \\ \text{Primary carbon}:& &\text{POA(50%)},& &\text{BC(50%)}. \nonumber \end{align} No initial SOA and $\rm SO_4$ are specified for the primary carbon mode because that mode only contains BC, POA and MOA. By specifying the $T$, $p$, $R$, $N_{_{t,i}}$, $D_{_{gn,d,i}}$, $\sigma_{_{g,i,}}$, $\rho_{_{l}}$, $M_{_{w,l}}$ and $f_{_{m,l,i}}$, we could then calculate the initial $q_{_{m,l,i}}$ as: \begin{equation} q_{_{m,l,i}} = \frac{\frac{\pi}{6} N_{_{t,i}} D_{_{gn,d,i}}^3 \exp \left( \frac{9}{2} \ln^2 \sigma_{_{g,i}} \right)} {\sum\limits_{l}\frac{f_{_{m,l,i}}}{\rho_{_{l}}}} \cdot f_{_{m,l,i}} \cdot \frac{\rho_{_{l}}RT}{1000 M_{_{w,l}} p} \,. \end{equation} \item Initial mass mixing ratio of $q_{_{v,SOAG}}$, which is set to $5 \times 10^{-10}$, close to the global mean value at the surface layer. \end{enumerate}

With these parameters given, totally 5 (= 4 (modes for aerosols + 1 (gas species)) verification tests are performed.

3.4 PASS/FAIL criterion

We first calculate the MAM solutions with different $\Delta t_{_{cond}}$ based on Eq. \eqref{sol-soag} and \eqref{sol-soa} (denoted as $q_{_{v,MAM}}$ and $q_{_{m,MAM,i}}$). We then calculate the reference solution with $\Delta t_{_{cond}} = 0.01$s (denoted as $q_{_{v,REF}}$ and $q_{_{m,REF,i}}$). Thus the relative error $\phi$ can be calcualted by:

\begin{align} \phi_{_v} &= \frac{\left| q_{_{v,MAM}} - q_{_{v,REF}} \right|}{q_{_{v,REF}}} \,, \\ \phi_{_{m,i}} &= \frac{\left| q_{_{m,MAM,i}} - q_{_{m,REF,i}} \right|}{q_{_{m,REF,i}}} \,. \end{align}

A self-convergence test is performed by plotting $\phi$ against $\Delta t_{_{cond}}$ on a log-log scale figure. Denoting $\tau = \frac{1}{\sum\limits_{i=1}^I C_{_{SOAG,i}}(t_{_0})}$ as the characteristic time scale of SOAG/SOA condensation, we further fit a line to the pairs of ($\Delta t_{_{cond}}$, $\phi$) when $\Delta t_{_{cond}} \le f_{_{\tau}}\tau$. $f_{_{\tau}}$ is user-specified number and it is determined based on the criterion that the solution starts to convege when $\Delta t_{_{cond}} \le f_{_{\tau}}\tau$.

An overall PASS is issued if the slope of fitted line is 1.0$\pm \delta$ for both $\phi_{_v}$ and $\phi_{_{m,i}}$ (Here $\delta$ is a user-specified parameter so that 1.0$\pm \delta$ is thought close enough to 1.0).

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/soag\_convergence}$ and the yaml file is under $\rm \color{blue}{tests/mam/restructed\_tests/yaml\_files\_for\_test\_runner/condensation/soag\_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 MAM and reference solutions for $q_{_{v,SOAG}}$ and $q_{_{m,SOA,i}}$. Assuming that the netcdf file is put under the "postprocess_input" folder, the following cells would calculate the relative errors $\phi$, perform a self-convergence test of MAM solutions and check whether a PASS/FAIL should be issued.

In [1]:
# load the library
import numpy as np
import xarray as xr
import postprocess_scripts.plot_conden_convergence as ps
In [2]:
# define the file path and name
nc_file   = "./postprocess_input/condensation_verification_test.nc"
outpath   = "./postprocess_output/"
In [3]:
# make the self-convergence plot of MAM solution for SOA
ps.plot_convergence(nc_file, outpath, var_type='aer', porder=[1], 
                    xlabel=np.arange(-2,5,1), ylabel=np.arange(-8,4,1), xlim=[-2,4], ylim=[-8,0],
                    left=0.1, right=1., top=1., bottom=0.1, hspace=0, wspace=0.1, 
                    nrows=1, ncols=1, fontheight=16)
In [4]:
# make the self-convergence plot of MAM solution for SOAG
ps.plot_convergence(nc_file, outpath, var_type='gas', porder=[1], 
                    xlabel=np.arange(-2,5,1), ylabel=np.arange(-6,2,1), xlim=[-2,4], ylim=[-6,0],
                    left=0.1, right=1., top=1., bottom=0.1, hspace=0, wspace=0.1, 
                    nrows=1, ncols=1, fontheight=16, loc='upper left')

If a user uses the same settings in Section 3.3 and does not change the source codes, the self-convergence results of SOAG/SOA would look like:

In [5]:
ps.show_figures("./trusted_output/",left=0.1,right=2.05,top=2.05,bottom=0.1,hspace=0.,wspace=0.)
In [6]:
# show the tau of SOAG/SOA condensation
tau = ps.calc_tau(nc_file)

τ of SOAG/SOA condensation: 1.9657e+03

Based on the plots of self-convergence test for SOAG/SOA above, both $\phi_{_v}$ and $\phi_{_{m,i}}$ start to converge when $\Delta t_{_{cond}} \le 100$s. Therefore, we specify $f_{_{\tau}} = \frac{1}{20}$ and fit a line to the pairs of ($\Delta t_{_{cond}}$, $\phi$) when $\Delta t_{_{cond}} \le f_{_{\tau}}\tau$. The resulted slope of fitted line is:

In [7]:
# show the slope of fitted line when dt <= ftau * tau
slope_gas, slope_aer = ps.calc_slope(nc_file, tau, ftau=1./20.)

slope of fitted line for SOAG:          1.01316
slope of fitted line for SOA in mode 1: 1.01317
slope of fitted line for SOA in mode 2: 1.01316
slope of fitted line for SOA in mode 3: 1.01318
slope of fitted line for SOA in mode 4: 1.01312

With the slopes known, we could then check whether a PASS/FAIL should be issued. Here we set $\delta = 0.02$ as a threshold (meaning the verification test passes if and only if the fitted slopes are all within 1$\pm$0.02 ).

In [8]:
# run this cell to issue PASS/FAIL
ps.issue_flag(nc_file, slope_gas, slope_aer, delta=0.02)

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] 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.

[3] 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.

[4] Daniel Zwillinger, ``CRC standard mathematical tables and formulae'', 2002.

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