Siyuan Wang | Research - Ch 9. SOA Formation

Ch 9. SOA Formation from Toluene

OVERVIEW

- In the previous chapter we learned how one can model condensational growth of aerosols using "moving bin" approach. But that approach was for non-volatile compounds (e.g. H2SO4) only. In this chapter we'll cover the condensation and evaporation of semi-volatile compounds!
- Depending on the volatility (described by saturation vapor pressure or saturation vapor concentration), a semi-volatile compound can exist partially in the particle-phase and partially in the gas-phase, and the condensation and evaporation together forms an equilibrium. Many models assume they'll reach equilibrium instaneously. However, numerious studies have shown that the equilibrium time scale can be long, i.e. it actually takes a while to reach equilibrium. In this chapter, the condensation and evaporation are treated in a kinetic approach which is computationally more expensive but more accurate.
- In the ambient environment, many primarily emitted volatile organic compounds (e.g. isoprene) are too volatile to form aerosols. However, these VOCs can undergo chemical reactions in the atmosphere, producing compounds with lower volatilities which can condense onto pre-existing aerosols. These are usually referred to as "secondary organic aerosols", or SOA. SOA plays a key role in the atmosphere. In this chapter, we'll use the model to explore the SOA formation from toluene, and we'll see how condensation and evaporation affect the aerosol size distribution.

1. Get Started

The prognostic equations for condensation and evaporation are fairly similar to what we had in the previous chapter:

dCoa/dt = kmt*(Cg - Cstar*Coa/total_Coa)

dCg/dt = -kmt*(Cg - Cstar*Coa/total_Coa)

where Cg and Coa are the gas-phase and particle-phase concentrations of a compound, kmt is again the phase-transfer coefficient - we'll use the same one from the previous Chapter, which involves Knudsen number. Cstar is the saturation concentration which can be calculate from saturation vapor pressure. Note that Kelvin effect can affect vapor pressure for small particles. If you have multiple aerosol bins, you need the first equation for each bin.

Go to the MCM website and generate the CRI mechanism file for toluene. I have the file uploaded in GitHub too (cri_export_toluene.eqn).
Go to GitHub and download the tutorial file (YAMCHA_tutorial_ch9). Drop the toluene KPP file in Google Colab!
Run the first 3 Cells. Load the KPP file in Cell #4. This mechanism has 259 reactions.

2. Set up the Seed Aerosols

Same as the previous chapter. We'll set an initial aerosol size distribution to start with. The initial size distribution is lognormal but you can use whatever shape you want. The initial distribution is shown on the right, as you can see we have 15 size bins ranging from 0.001 um to 10 um.

N0 = 7e+5 # total aerosol number concentration, cm^-3

Dp_mean_um = 0.012 # initial median diameter

sigma = 2.7 # aerosol size dist sigma

N_bins = 15

init_Dp_bounds_um = 10**np.linspace(-3, 1, N_bins+1)

Then proceed to Cell #6. Here we'll set SOA formation related chemical reactions.

TOLUENE + OH = TOLUENE + OH + 0.1364*SOAG0 + 0.0101*SOAG1 + 0.0763*SOAG2 + 0.2157*SOAG3 + 0.0738*SOAG4

This is the main reaction for toluene+OH, which produces 5 semi-volatile species and they all can partition into the aerosol-phase (depending on their volatilities): SOAG0-SOAG4, where "G" means they're in the gas-phase. They have different molar yields (SOAG3 has the highest molar yield).

This is the so-called "volatility basis set (VBS)" approach, which lumps all the semi-volatile species into different volatility bins based on their individual volatilities. As a result, SOAG0-SOAG4 have saturation vapor concentrations of 0.01, 0.1, 1, 10, 100 ug/m3 (at 300K). See Tilmes et al. 2019.
Note here Toluene and OH are regenerated. This is because the main mechanism (the KPP file you loaded, cri_export_toluene.eqn), the toluene chemistry (including it's initial OH attack) is taken care of, so if you don't regenerate toluene and OH in this reaction, you'll be double-counting the toluene and OH removal.
We'll also have a few reactions taking care of the wall loss of these species. At this point let's assume they're irreversibly removed with a first order lost rate coefficient (kwall).
Don't worry about all the "None". This is because when the chemical mechanism is loaded from MCM-KPP files, a RO2 list is attached to the mechanism. When attaching additional reactions to the mechanism, the RO2 list will be extended but with no additional elements.

3. Set up the Condensation/Evaporation

This is a main module where we'll "fake" the kinetic condensation and evaporation as chemical reactions so they can be solved using the ODE solver. The toluene oxidation produces 5 semi-volatile species, SOAG0-4, and each can partition into all 15 aerosol size bins. That is:

SOAG[N] = SOA[N,i]

SOA[N,i] = SOAG[N]

where N is the index of the semi-volatile species in the gas-phase, i.e. N = 0-4 (remember we have 5 of them, SOAG0-SOAG4). i represents the i-th aerosol bin. So SOA[N,i] is the corresponding aerosol-phase species of SOAG[N] in the i-th aerosol bin. I have a code snippet that'll generate all these "reactions":

RN = 0

tmp_rindexlist,tmp_reactionlist,tmp_coefflist,tmp_ro2 = [],[],[],[]

for iff in [0,1,2,3,4]:

ffgspc = 'SOAG' + str(iff)

ffaspc = 'SOA' + str(iff)

for ibin in range(len(init_bin_center_um)):

for itmp in [0,1,2]:

# --- condensation: SOAGn = SOAn_i (i: ith aerosol bin)

if itmp==0:

tmp_rindexlist.append('P'+str(RN))

tmp_reactionlist.append(ffgspc+' = '+ffaspc+'_bin'+str(ibin))

tmp_coefflist.append('kmt[%d]' % (ibin))

tmp_ro2.append(None)

RN += 1

# --- evaporation: SOAn_i = SOAGn (i: ith aerosol bin)

if itmp==1:

tmp_rindexlist.append('P'+str(RN))

tmp_reactionlist.append(ffaspc+'_bin'+str(ibin)+' = '+ffgspc)

tmp_coefflist.append('kmt_Cstar_f[%d,%d]' % (iff,ibin))

tmp_ro2.append(None)

RN += 1

ChemMech_GasParticle = pd.DataFrame(zip(tmp_reactionlist,tmp_coefflist,tmp_ro2))

ChemMech_GasParticle.index = tmp_rindexlist

ChemMech_GasParticle.columns = ['reaction', 'rate_coefficient','RO2']

ChemMech_GasParticle = ChemMech_GasParticle.T

ChemMech_GasParticle.T

Then we'll have anothe line that merges all of these processes into one big mechanism which will be loaded into the model.

# ==========================

# Merge into a big mechanism

# ==========================

ChemMech_merge = pd.concat([ChemMech,ChemMech_SOA,ChemMech_GasParticle],axis=1)

ChemMech_merge.T

If all goes right, a summary of the merged mechanism will be printed, which has 415 rows, i.e. 415 reactions.

4. Main Model Cell

Proceed to Cell #8 which is the main model cell. Firstly we need to set the phase-transfer coefficient that is used to calculate the gas-particle transfer, which is the same function that we used in the previous chapter:

def kmt_gas_aerosol(in_radius_cm, in_Dg_cm2_s, in_N_cm3, in_alpha,

in_temp, in_press, in_MW_gas_kg_mol):

MW_air_kg_mol = 0.029

# DynamicViscosity_kg_m_s = 1.8325e-5*((416.16/(in_temp+120.))*(in_temp/296.16)**1.5)

DynamicViscosity_kg_m_s = 1.715747771E-5 + 4.722402075E-8 * (in_temp-273.15) \

- 3.663027156E-10 * ((in_temp-273.15)**2.0) \

+ 1.873236686E-12 * ((in_temp-273.15)**3.0) \

- 8.050218737E-14 * ((in_temp-273.15)**4.0)

thermalspeed_m_s = np.sqrt(8.0*8.314*in_temp/np.pi/in_MW_gas_kg_mol)

rho_air_kg_m3 = in_press*MW_air_kg_mol/8.314/in_temp

MeanFreePath_m = 3.*DynamicViscosity_kg_m_s/rho_air_kg_m3/thermalspeed_m_s

Kn = (MeanFreePath_m*100.) / in_radius_cm

f = 0.75*in_alpha*(1.+Kn)/( Kn*(1.+Kn) + 0.283*in_alpha*Kn + 0.75*in_alpha )

return 4.*np.pi*in_radius_cm*in_Dg_cm2_s*in_N_cm3*f

Most of the stuff in this cell should look very familiar to you. We'll set the total run time to 18 hours, and we'll output the results every 600 seconds (10 min). Make sure that the merged mechanism (ChemMech_merge) is loaded! A few important things:
Photolysis: we discussed this in Chapter 3. Here we'll use photolysis frequencies parameterized as a function of solar zenith angle (SZA), which is a good approximation for ambient conditions but not necessarily for chamber conditions. Here we'll just tweak SZA to mimic the photolysis in the chamber condition.
Initial condition: here we'll have low RH and low NOx, so H2O2 photolysis will be the main radical source. The initial concentrations of toluene and H2O2 are set to 37.6 ppb and 2 ppm, respectively.
Wall loss: this is important. Here we'll set up the vapor wall loss of those semi-volatile species produced from toluene (SOAG0-4), and their first order wall loss rate coefficient (kwall) is set to 2.5e-4 s^-1. Technically you also should setup wall loss of aerosols (e.g. seed aerosols and SOA produced) but for simplicity let's ignore those for now.
Volatilities of the semi-volatile species: because we're using VBS, these are set values. But if you want to deal with other species (e.g. condensation/evaporation of water), you need to get their saturation vapor pressure or saturation vapor concentrations.
Kelvin effect: this affects the vapor pressure over small particles. We'll use the surface tension for water which supposedly works also for many applications. Kelvin effect depends on aerosol size, so you have to recalcualte the Kelvin effect correction factor at each time step when the aerosol sizes are updated.
Gas-particle partition kinetics: this is updated also at the end of each chemistry time step.

# --- update size and kinetics

for ibin in range(len(init_bin_center_um)):

# --- calculate total SOA in each size bin, also total SOA per particle

tot_OA_per_particle_g, tot_OA_per_bin_ug_m3 = 0., 0.

for ivbs in range(len(Cstar_volatility_bins_ug_m3)):

tot_OA_per_particle_g += conc[i,species_list.index('SOA'+str(ivbs)+'_bin'+str(ibin))]/init_bin_N_cm3[ibin]/6.0232e+23*MW_g_mol

tot_OA_per_bin_ug_m3 += conc[i,species_list.index('SOA'+str(ivbs)+'_bin'+str(ibin))]*MW_g_mol/6.0232e+11

# --- new aerosol mass: per particle. seed mass + SOA mass

m_per_particle_g[ibin] = seed_ug_m3[ibin]*1e-6/1000000./init_bin_N_cm3[ibin] + tot_OA_per_particle_g

# --- get new aerosol size

new_Dp_um[ibin] = 2.*1e+4*((m_per_particle_g[ibin]/rho_g_cm3)*(3./4.)/np.pi)**(1./3.)

# --- Kelvin effect

if if_Kelvin: KelvinFactor = np.exp(4.*MW_g_mol*surface_tension_water_N_m/8.314/temp/(rho_g_cm3*1000000.)/(new_Dp_um[ibin]*1e-6))

else: KelvinFactor = 1.

# --- update kinetics

kmt[ibin] = kmt_gas_aerosol(0.5*new_Dp_um[ibin]*1e-4, 0.1, init_bin_N_cm3[ibin], mass_accommodation_coeff, temp, press, MW_g_mol/1000.)

for ivbs in range(len(Cstar_volatility_bins_ug_m3)):

kmt_Cstar_f[ivbs,ibin] = kmt[ibin] * KelvinFactor*Cstar_volatility_bins_ug_m3[ivbs]/MW_g_mol / (seed_ug_m3[ibin]/seed_MW_g_mol+tot_OA_per_bin_ug_m3/MW_g_mol)

As you can see, it loops over each aerosol bin, and for each bin we firstly calculate the total SOA mass per particle which is added to the seed mass to update the size (diameter). Then with the new aerosol size, we'll calculate the Kelvin effect correction factor. Lastly, the kinetics (kmt) is updated.

5. Run the Model!

Let's run the model! In the end a few plots are generated, and the results are attached below.
The left plot shows the decay of toluene and the formation of the 5 semi-volatile species (SOAG0-4). We can see toluene decays to ~10 ppb at the end of the 18 hour simulation. The concentrations of SOAG0-4 remain low, because they're both partitioned into the particle-phase, also deposited onto the walls. The middle plot shows the OA formation, with total OA showing as a thick black line, reaching ~45 ug/m3 at the end of the simulation. The colorful lines show each of the semi-volatile component: the blue line (SOA4) is one of the highest, because SOAG0 has the lowest volatility (C*=0.01 ug/m3) and the second highest molar yield. The red line (SOA3) is one of the highest, although it's somewhat volatile (C* = 10 ug/m3) but it has the highest molar yield. The right plot shows the gas-phase fraction of each component. We can see that the most volatile species, SOAG4 has the highest gas-phase fraction, while SOAG0 and SOAG1 have the lowest gas-phase fractions because of their low volatilities. It's also important to note from the right plot that the partitioning of all species do not reach equilibrium after 18-hour simulation.

Since we're kinetically modeling the evoluation of the aerosol size distribution, we can plot the size distributions and see how they change! The plot on the right shows the dN/dlogDp and dV/dlogDp at the beginning and the end of the simulation. We can see that the number distribution becomes larger and narrower because of the condensational growth. The change in volume distribution also grows taller but the shape doesn't change as much. Remember in this simulation we ignore the wall loss of aerosols, and we don't have coagulation as well.
Also shown here is the "banana" plot, which is the number distribution (dN/dlogDp) evolution as a function of time. The particle growth is very interesting!

6. Kelvin Effect

Kelvin effect describes the vapor pressure over a curved surface (e.g. an aerosol particle) always exceeds the vapor pressure over a flat surface. The top figure on the right demonstrates the Kelvin effect for water. As you can see, on a curved surface, there are "fewer" molecules adjacent, so it's easier for the molecule to "escape".
The middle and bottom figures show the aerosol number size distribution with Kelvin effect turned on and off (in the model, Kelvin effect is controlled by a True/False switch). You can see that if Kelvin effect is turned on, the growth of very small particles (e.g. a few to 10 nm) is greatly limited. Note the different y-axis scales on these two plots.
In the atmosphere, Kelvin effect has major impact on nucleation. For aerosols larger than a few hundreds of nm, the Kelvin effect is small. For cloud/fog droplets, the Kelvin effect could probably be ignored.

7. Vapor Wall Loss

It is known that in chamber experiments, aerosol particles can lost to chamber walls. The semi-volatile vapors (e.g. SOAG0-4 in this model) can also deposit on chamber walls, representing an important sink for these vapors. This was not known until somewhat recently (Zhang et al. 2014). If Vapor wall loss is not considered, the SOA yields derived from chamber studies will be biased.
Here we can run a set of simulation with vapor wall loss turned off (set kwall to zero) and examine how it'll affect the results. The plot is attached on the right. As you can see, if vapor wall loss is turned off, the simulated total organic aerosol production is increased by ~3x.
Studies have shown that the vapor wall loss is actually an irreversible process. For simplicities, in this demo the wall loss is treated as an irreversible process but there are methods out there if you want to properly model the reversible wall loss.

8. Summary

In this chapter, we learned some fundamentals of the condensation and evaporation of semi-volatile compounds, which in principle is similar to what we had in the previous chapter. Method-wise, the treatment of the condensation and evaporation is very similar to the gas-aqueous partitioning that we showed in Chapter 7!
We also coded the VBS scheme for toluene and solve that using the kinetic approach rather than the equilibrium approach. The VBS framework is widely used in research. This chapter only covered SOA formation from toluene. Similarly, SOA can be produced from a wide variety of other VOCs, e.g. isoprene, pinenes. The VBS framework can be used for other VOC precursors.
This model is loosely configured based on the chamber study reported in Zhang et al. 2014. As you can tell, I'm not trying very hard to mimic the exact conditions reported in that study but if you want you can certainly refine the model!

Google Sites

Report abuse