Aerosol condensation process and its verification test: correct implementation of semi-analytical solutions of $\rm H_2SO_4/SO_4$¶
$\color{red}{\text{Last update: 09-09-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 $\rm H_2SO_4$ 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} \\ \frac{dq_{_{v,H_2SO_4}}}{dt} &= P_{_{chem}} - q_{_{H_2SO_4}} \sum_{i=1}^{I} C_{_{H_2SO_4,i}} \,, \label{qgh2so4} \\ \frac{dq_{_{m,SO_4,i}}}{dt} &= C_{_{H_2SO_4,i}} \times q_{_{v,H_2SO_4}} \,, i = 1,2,...,I \,, \label{qaso4} \end{align}where the definitions of variables in Eq. \eqref{kmt-h2so4} can be found in a separate model documentation; $q_{_{v,H_2SO_4}}$ is the mass mixing ratio of $\rm H_2SO_4$ (unit: kmol of $\rm H_2SO_4$ per kmol of dry air); $P_{_{chem}}$ is the net chemical production rate of $\rm H_2SO_4$ due to the gas- and aqueous-phase chemistry (unit: kmol of $\rm H_2SO_4$ per kmol of dry air per second); $q_{_{m,SO_4,i}}$ is the mass mixing ratio of $\rm SO_4$ in mode $i$.
To solve for $q_{_{v,H_2SO_4}}$ and $q_{_{m,SO_4,i}}$, MAM assumes that $P_{_{chem}}$ and $C_{_{H_2SO_4,i}}$ are fixed within one time step $\Delta t_{_{cond}}$ ($\Delta t_{_{cond}} = t_{_1} - t_{_0}$). Therefore, semi-analytical solutions can be derived for $q_{_{v,H_2SO_4}}$ and $q_{_{m,SO_4,i}}$. More details about the basics notations, model parameterization, assumptions and numerical implementation can be found in a separate ''MAM_Basics'' and model documentation.
2 Goals¶
Though there are semi-analytical solutions for $q_{_{v,H_2SO_4}}$ and $q_{_{m,SO_4,i}}$, we do not know whether the numerical implementation is correct (e.g., no bug in the index of array or array name). Thus in this test, we aim at examing the correct implementation of semi-analytical solutions of $\rm H_2SO_4/SO_4$ in MAM. Note that in this test, $P_{_{chem}}$ and $C_{_{H_2SO_4,i}}$ are fixed over the whole simulation time and we do not consider the existence of SOAG/SOA.
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¶
In this test, we use typical number concentration of aerosols in the atmosphere and $\Delta t_{_{cond}} = 1800s$ as default in MAM. Therefore, MAM solution takes the semi-analytical formulae for $\rm H_2SO_4/SO_4$, based on the assumptions mentioned in Section 1. The semi-analytical solutions of Eq. \eqref{qgh2so4} and \eqref{qaso4} at $t = t_{_1}$ read:
\begin{align} q_{_{v,H_2SO_4}}(t_{_1}) &= \left( q_{_{v,H_2SO_4}}(t_{_0}) - \frac{P_{_{chem}} (t_{_0})} {C_{_{sum,H_2SO_4}} (t_{_0})} \right) \exp^{-C_{_{sum,H_2SO_4}} (t_{_0}) \Delta t_{_{cond}}} + \frac{P_{_{chem}} (t_{_0})} {C_{_{sum,H_2SO_4}} (t_{_0})} \,, \label{semi-ana-h2so4} \\ q_{_{m,SO_4,i}}(t_{_1}) &= q_{_{m,SO_4,i}}(t_{_0}) + \frac{C_{_{H_2SO_4,i}} (t_{_0})} {C_{_{sum,H_2SO_4}} (t_{_0})} \left[ \left( q_{_{v,H_2SO_4}}(t_{_0}) - \frac{P_{_{chem}} (t_{_0})} {C_{_{sum,H_2SO_4}} (t_{_0})} \right) \left( 1-\exp^{-C_{_{sum,H_2SO_4}} (t_{_0}) \Delta t_{_{cond}}} \right) + P_{_{chem}} (t_{_0}) \Delta t_{_{cond}} \right] \,. \label{semi-ana-so4} \end{align}where $C_{_{sum,H_2SO_4}} (t_{_0}) = \sum\limits_{i=1}^{I}C_{_{H_2SO_4,i}} (t_{_0})$.
3.2 Reference solution¶
Define the following quantity $q_{_{tot}}(t)$:
\begin{equation} \label{qtot} q_{_{tot}}(t) = q_{_{v,H_2SO_4}}(t_{_0}) + \sum_{i=1}^I q_{_{m,SO_4,i}}(t_{_0}) + P_{_{chem}} (t-t_{_0}) \,. \end{equation}We could rewrite the governing equations of Eq. \eqref{qgh2so4} and \eqref{qaso4} as:
\begin{align} \frac{dq_{_{m,SO_4,i}}}{dt} &= C_{_{H_2SO_4,i}} \left( q_{_{tot}}(t) - \sum_{i=1}^I q_{_{m,SO_4,i}} \right) = F_{_i}(t, q_{_{m,SO_4,1}},..., q_{_{m,SO_4,I}}) \,, i=1,2,...I, \label{qaso4-rk4} \\ q_{_{v,H_2SO_4}} &= q_{_{tot}}(t) - \sum_{i=1}^I q_{_{m,SO_4,i}} \,. \label{qgh2so4-rk4} \end{align}If we define the vectors $\mathbf{q}_{_{m,SO_4}} = (q_{_{m,SO_4,1}},..., q_{_{m,SO_4,I}})^T$ and $\mathbf{F} = (F_{_1},..., F_{_I})^T$, the Eq. \eqref{qaso4-rk4} can be further rewritten as:
\begin{equation} \label{qaso4-rk4v2} \frac{d\mathbf{q}_{_{m,SO_4}}}{dt} = \mathbf{F}(t, \mathbf{q}_{_{m,SO_4}}) \end{equation}We apply the Runge-Kutta Fourth-order (RK4) method to solve Eq. \eqref{qaso4-rk4v2} for $\mathbf{q}_{_{m,SO_4}}$. The solution of $\mathbf{q}_{_{m,SO_4}}(t)$ at $t = t_{_1}$ reads:
\begin{align} \mathbf{K}_{_1} &= \mathbf{F}(t, \mathbf{q}_{_{m,SO_4}}) \,, \\ \mathbf{K}_{_2} &= \mathbf{F}(t+\frac{1}{2}\Delta t_{_{cond}}, \mathbf{q}_{_{m,SO_4}}+\frac{1}{2}\Delta t_{_{cond}} \mathbf{K}_{_1}) \,, \\ \mathbf{K}_{_3} &= \mathbf{F}(t+\frac{1}{2}\Delta t_{_{cond}}, \mathbf{q}_{_{m,SO_4}}+\frac{1}{2}\Delta t_{_{cond}} \mathbf{K}_{_2}) \,, \\ \mathbf{K}_{_4} &= \mathbf{F}(t+\Delta t_{_{cond}}, \mathbf{q}_{_{m,SO_4}}+\Delta t_{_{cond}} \mathbf{K}_{_3}) \,, \\ \mathbf{q}_{_{m,SO_4}}(t_{_1}) &= \mathbf{q}_{_{m,SO_4}}(t_{_0}) + \frac{\Delta t_{_{cond}}}{6} \left( \mathbf{K}_{_1} + 2\mathbf{K}_{_2} + 2\mathbf{K}_{_3} + \mathbf{K}_{_4} \right) \,. \end{align}After solving $\mathbf{q}_{_{m,SO_4}}(t_{_1})$, $q_{_{v,H_2SO_4}}(t_{_1})$ can be calculated by Eq. \eqref{qgh2so4-rk4}.
In this verification test, we would generate a series of solutions with different $\Delta t_{_{cond}}$. A self-convergence test is performed to show the confidence of solutions solved by RK4 method (expected $\rm 4^{th}$-order convergence).
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} \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}}$ for MAM and reference solutions: \begin{itemize} \item $\Delta t_{_{cond}} = 1800$ is used to generate the MAM solution; \item $\Delta t_{_{cond}} \in $ [1800.,900.,600.,300.,100.,60.,30.,15.,5.,1.,0.1] are used to perform a self-convergence test of the solutions solved by the RK4 method and determine the reference solutions. \end{itemize} \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{$\rm SO_4$(50%)}, \nonumber \\ \text{ Aitken}:& &\text{SS(50%)},& &\text{$\rm SO_4$(50%)}, \nonumber \\ \text{ Coarse}:& &\text{DST(50%)},& &\text{$\rm SO_4$(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,H_2SO_4}}$, which is set to $10^{-13}$, close to the global mean value at the surface layer. \item The net chemical production rate $P_{_{chem}}$ for $\rm H_2SO_4$, which is set to $10^{-16}$ as an example. \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 use the MAM solution at $\Delta t = 1800$s as the reference solution and check the convergence of RK4 solutions (for both $\rm H_2SO_4$ and $\rm SO_4$) from $\Delta t = 0.1$s to $\Delta t = 1800$s. An overall PASS is issued if (i) a $4^{th}$-order convergence is observed in a log-scale plot of $\Delta t$ against relative error within a certain range of $\Delta t$ (say [$\rm dt_{lo}$, $\rm dt_{hi}$]) and (ii) the relative error does not reduce or even increases when $\Delta t$ becomes small enough. In this test, a $4^{th}$-order convergence is reached when the slopes of fitted lines over [$\rm dt_{lo}$, $\rm dt_{hi}$] are within $4\pm\delta$, where $\delta$ is a user-specified constant.
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/h2so4\_sol\_correctness}$ and the yaml file is under $\rm \color{blue}{tests/mam/restructed\_tests/yaml\_files\_for\_test\_runner/condensation/h2so4\_sol\_correctness.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,H_2SO_4}}$ and $q_{_{m,SO_4,i}}$. Assuming that the netcdf file is put under the "postprocess_input" folder, the following cells would perform a self-convergence test of reference solutions, calculate the relative errors, and check whether a PASS/FAIL should be issued.
5.1 Confidence of RK4 solutions¶
# load the library
import numpy as np
import xarray as xr
import postprocess_scripts.plot_conden_convergence as ps
import sys
import warnings
# turn off warning message
if not sys.warnoptions:
warnings.simplefilter("ignore")
# define the file path and name
nc_file = "./postprocess_input/condensation_verification_test.nc"
outpath = "./postprocess_output/"
# make the convergence plot of RK4 solution for SO4
ps.plot_convergence(nc_file, outpath, var_type='aer', porder=[4],
xlabel=np.arange(-1,5,1), ylabel=np.arange(-16,4,2), xlim=[-1,4], ylim=[-16,2],
left=0.1, right=1., top=1., bottom=0.1, hspace=0, wspace=0.1,
nrows=1, ncols=1, fontheight=16)
# make the convergence plot of RK4 solution for H2SO4
ps.plot_convergence(nc_file, outpath, var_type='gas', porder=[4],
xlabel=np.arange(-1,5,1), ylabel=np.arange(-14,2,2), xlim=[-1,4], ylim=[-12,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')
Based on the results above, we obtain a $4^{th}$-order convergence rate for both $\rm H_2SO_4$ and $\rm SO_4$, which indicates the correct implementation of RK4 method.
5.2 Use MAM solution as refence solution and re-check the convergence of RK4 solutions¶
After we verify the correct implementation of RK4 method, we use the MAM solution as the reference solution and re-check the convergence of RK4 solutions. Here we include the RK4 solutions from $\Delta t = 0.1$s to $\Delta t = 1800$s. If the implementation of MAM solution is correct, we should again observe a $4^{th}$-order convergence rate as shown in Section 5.1.
# make the convergence plot of RK4 solution for SO4, using MAM solution as reference
ps.plot_convergence(nc_file, outpath, var_type='aer', porder=[4], use_mam='TRUE',
xlabel=np.arange(-2,5,1), ylabel=np.arange(-16,4,2), xlim=[-2,4], ylim=[-16,2],
left=0.1, right=1., top=1., bottom=0.1, hspace=0, wspace=0.1,
nrows=1, ncols=1, fontheight=16)
# make the convergence plot of RK4 solution for H2SO4, using MAM solution as reference
ps.plot_convergence(nc_file, outpath, var_type='gas', porder=[4], use_mam='TRUE',
xlabel=np.arange(-2,5,1), ylabel=np.arange(-14,2,2), xlim=[-2,4], ylim=[-14,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')
5.3 Check the convergence and issue PASS/FAIL¶
We check the convergence rate (the slope of fitted line using the data points between $\Delta t = 30$s and $\Delta t = 600$s in Section 5.2).
# calculate the slopes of fitted lines
slope_gas, slope_aer = ps.check_convergence(nc_file, 600., 30.)
Set $\delta$ as 0.1, and we could issue PASS/FAIL by running the following cell:
# run this cell to issue PASS/FAIL
delta = 0.1
if np.abs(slope_gas - 4) < delta and np.all(np.abs(slope_aer - 4) < delta):
print("PASS")
else:
print("FAIL")
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] , ``Python function to generate the Gauss–Hermite quadrature points'', . online