Aerosol condensation process and its verification test: self-convergence rate of semi-implicit solutions of SOAG/SOA against different $\alpha_{_{astem}}$

$\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

MAM implemented a dynamical substepping approach within one $\Delta t_{_{cond}}$ to solve for the semi-implicit solutions of SOAG/SOA. The varying-length substep $\Delta t_{_{SOAG}}$ is controlled by a user-specified error-control parameter $\alpha_{_{astem}}$ with smaller $\alpha_{_{astem}}$ leading to smaller $\Delta t_{_{SOAG}}$. In this test, we aim at examining whether the relative errors of solutions decrease linearly with the decrease of $\alpha_{_{astem}}$ in a log-log scale plot (e.g., $1^{st}$-order convergence rate). We would perform a self-convergence test of semi-implicit solutions of SOAG/SOA against different $\alpha_{_{astem}}$ 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$ at this moment.

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}

The varying-length substep $\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}$) is first calculated by:

\begin{align} q_{_{v,SOAG,sat,i}} (t_{_{old}}) &= \frac{q_{_{m,SOA,i}}(t_{_{old}})}{\max \left( q_{_{m,OPOA,i}}(t_{_0}) + q_{_{m,SOA,i}}(t_{_{old}}), 10^{-20} \right)} q_{_{v,SOAG,equ}}(t_{_0}) \,, \\ \Phi_{_i} (t_{_{old}}) &= \frac{q_{_{v,SOAG}}(t_{_{old}}) - q_{_{v,SOAG,sat,i}}(t_{_{old}})}{\max \left( q_{_{v,SOAG}}(t_{_{old}}), q_{_{v,SOAG,sat,i}}(t_{_{old}}), 10^{-20} \right)} \,, \\ \Delta t_{_{SOAG}} &= \min \left( \frac{\alpha_{_{astem}}}{\sum\limits_{i=1}^{I} \left( C_{_{SOAG,i}}(t_{_0}) \left| \Phi_{_i}(t_{_{old}}) \right| \right)}, t_{_1} - t_{_{old}} \right) \,. \label{dtsub} \end{align}

The semi-implicit solutions of Eq. \eqref{qgsoag} and \eqref{qasoa} from $t_{_{old}}$ to $t_{_{new}}$ then 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.

In this test, the MAM solutions of SOAG/SOA are generated using the same $\Delta t_{_{cond}}$ but different $\alpha_{_{astem}}$ ranging from $10^{-4}$ to 1. The default value of $\alpha_{_{astem}}$ is 0.05 in MAM.

3.2 Reference solution

The reference solution is generated by the same method mentioned in Section 3.1 but with $\alpha_{_{astem}} = 10^{-5}$.

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}}$, which is set to 1800 seconds. \item Error-control paramerter $\alpha_{_{astem}}$ (unitless). In this test, we use the following values to generate the MAM solutions: \begin{equation} \alpha_{_{astem}} \in [1.,0.1,0.05,0.01,10^{-3},10^{-4}] \nonumber \end{equation} We use $\alpha_{_{astem}} = 10^{-5}$ 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 $\alpha_{_{astem}}$ 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 $\alpha_{_{astem}} = 10^{-5}$ (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 $\alpha_{_{astem}}$ on a log-log scale figure. We further fit a line to the pairs of ($\alpha_{_{astem}}$, $\phi$) where $\alpha_{_{astem}} \le 0.05$.

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\_substep\_convergence}$ and the yaml file is under $\rm \color{blue}{tests/mam/restructed\_tests/yaml\_files\_for\_test\_runner/condensation/soag\_substep\_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. Note that the asterisk symbol in the figures below means the relative error using the default $\alpha_{_{astem}} = 0.05$ in MAM.

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(-5,2,1), ylabel=np.arange(-8,4,1), xlim=[-5,1], ylim=[-8,0],
                    left=0.1, right=1.1, top=0.9, bottom=0.1, hspace=0, wspace=0.1, 
                    nrows=1, ncols=1, fontheight=20)
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(-5,2,1), ylabel=np.arange(-6,2,1), xlim=[-5,1], ylim=[-5,0],
                    left=0.1, right=1.1, top=0.9, bottom=0.1, hspace=0, wspace=0.1, 
                    nrows=1, ncols=1, fontheight=20, 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]:
# trusted output for SOAG/SOA 
ps.show_figures('./trusted_output/', left=0.1, right=2.2, top=2.2, bottom=0.1, hspace=0, wspace=0)

Based on the plots of self-convergence test for SOAG/SOA above, we further fit a line to the pairs of ($\alpha_{_{astem}}$, $\phi$) where $\alpha_{_{astem}} \le 0.05$. The resulted slope of fitted line is:

In [6]:
# show the slope of fitted line
slope_gas, slope_aer = ps.calc_slope(nc_file)

slope of fitted line for SOAG:          1.01121
slope of fitted line for SOA in mode 1: 1.01123
slope of fitted line for SOA in mode 2: 1.01121
slope of fitted line for SOA in mode 3: 1.01124
slope of fitted line for SOA in mode 4: 1.01116

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 [7]:
# 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