In previous chapters, we have learnt how to load KPP mechanism, as well as how to manually type-in mechanism. We also played around with VOC-NOx chemistry. Moreover, we set up the model to mimic the ambient conditions, such set diurnal variations for photolysis, mixing layer height, and depositions and emissions, etc.
In this chapter, we will briefly introduce heterogeneous chemistry and its treatment in the model. We will do this by playing around with a very interesting chemistry system: N2O5-ClNO2.
We will manually add a few more reactions. Burkholder et al JPL compication 2015 is a very good source for kinetics!!
Finally, we will take a look at some unique VOC signatures due to Cl oxidation!
N2O5 is an important nighttime NOx reservoir. N2O5 is also quite "sticky", i.e. when N2O5 hits a surface, the chance that N2O5 is "permanently" lost at the surface is quite high. Not all molecules striking a surface would trigger a chemical reaction. The probability of a molecule stick to a surface and triggers a chemical reaction is called reactive uptake coefficient.
So, when N2O5 hits a surface, what does it do? Well, if the N2O5 hits a "wet" surface (almost all surface has at least a teeny tiny amount of surfacial water) it would produce two HNO3 (which would usually rapidly disassociate into hydrogen ions and nitrate ions). But, if there is chloride on the surface, N2O5 may react with chloride and release ClNO2 into the gas-phase. In the ambient air, you'll amost always find at least some chloride.
ClNO2 is quite stable at night. In the morning, ClNO2 would undergo photolysis and produce Cl atoms, which is super reactive and can oxidize a lots VOCs, much faster than OH radicals. More interestingly, Cl oxidation has unique "signatures" that can be distinguished from OH oxidation.
Firstly go to the MCM website, and download the KPP mechanism containing these VOCs: Methane (CH4), ethane (C2H6), propane (C3H8), n-butane (n-C4H10), and i-butane (i-C4H10). Rename the file to "mcm_export_C1_C4.eqn". I also have this file in GitHub (link).
Grab the all-in-one Jupyter Notebook (YAMCHA_tutorial_ch5) from GitHub. Dump the KPP file in Colab.
Then run the first 3 Cells. In Cell #4, load the KPP file. You may need to edit the path to the file in order to load that correctly. If loaded correctly, a summary table will be printed, showing that the mechanism has 765 rows.
Then proceed to Cell #6, run it. This defines the sine-shaped function that we'll use later.
Now proceed to Cell #7, create a few additional reactions as well as a user-defined rate coefficient. Run Cells #7 and #8.
def k_termolecular(k0_300, n, kinf_300, m, temp, c_M):
k0_T = k0_300 * (temp/300.)**(-1.*n)
kinf_T = kinf_300 * (temp/300.)**(-1.*m)
return (k0_T*c_M/(1. + k0_T*c_M/kinf_T)) * 0.6**((1.+(np.log10(k0_T*c_M/kinf_T))**2.)**-1.)
ChemMech_add = pd.DataFrame({'D0': ['O3 = ', 'dep_rate_O3_sec', None],
'D1': ['NO2 = ', 'dep_rate_NO2_sec', None],
'D2': ['HNO3 = ', 'dep_rate_HNO3_sec', None],
'D3': ['N2O5 = ', 'dep_rate_N2O5_sec', None],
'H0': ['N2O5 = CLNO2', 'k_het_N2O5_CLNO2', None],
'H1': ['CLNO2 = CL + NO2', '0.005*jval[4]', None],
'H2': ['CL + NO2 = CLNO2', 'k_termolecular(1.3e-30, 2, 1e-10, 1, temp, M)', None],
})
ChemMech_add.index = ['reaction', 'rate_coefficient', 'RO2']
ChemMech_add.T
As you can see, here we add 4 *fake reactions* for the dry deposition for O3, NO2, HNO3, and N2O5. Then we add another 3 reactions:
N2O5 = CLNO2: this is a heterogeneous reaction happening on aerosol surfaces. We'll take care of its rate coefficient (k_het_N2O5_CLNO2) later.
CLNO2 = CL + NO2: this is the photolysis reaction of CLNO2. Its photlysis frequency is scaled to 0.5% of j-NO2 (jval[4]). You know how to estimate this? Simply integrate the absorption cross section of ClNO2 between ~290 and ~700 nm, and normalized to the absorption cross section of NO2.
CL + NO2 = CLNO2: this is one production pathway of CLNO2, and it's rate coefficient is given by a termolecular fomulation, defined above.
IMPORTANT: Python is case sensitive. In MCM, they do have CL radicals, and they use all uppercase letters for it, i.e. CL instead of Cl or cl.
As a result, in this model, use uppercase letters for CL!
Move further down to Cell #9:
# --- initial condition. need to convert to molec/cm3
conc_init['O3'] = (20.) *press*(6.0232E+8)/(8.314*temp)
conc_init['NO'] = (5.) *press*(6.0232E+8)/(8.314*temp)
conc_init['NO2'] = (10.) *press*(6.0232E+8)/(8.314*temp)
conc_init['CO'] = (100.) *press*(6.0232E+8)/(8.314*temp)
conc_init['CH4'] = (1800.) *press*(6.0232E+8)/(8.314*temp)
conc_init['C2H6'] = (5.) *press*(6.0232E+8)/(8.314*temp)
conc_init['C3H8'] = (2.) *press*(6.0232E+8)/(8.314*temp)
conc_init['NC4H10'] = (1.) *press*(6.0232E+8)/(8.314*temp)
conc_init['IC4H10'] = (1.) *press*(6.0232E+8)/(8.314*temp)
# --- set emissions and depositions
dep_vel_O3_cm_s = 0.2
dep_vel_NO2_cm_s = 0.2
dep_vel_HNO3_cm_s = 1.
dep_vel_N2O5_cm_s = 1.
Here we set the initial concentrations of a number of compounds, also the deposition velocities of O3, NO2, HNO3, and N2O5.
Move further down and you'll find this chunk of codes, which set the heterogeneous chemistry of N2O5.
# --- set N2O5-ClNO2 heterogeneous chemistry
Yield_ClNO2 = 0.6
gamma_N2O5 = 0.03
AerosolSurfaceArea_um2_cm3 = 200.
ThermalSpeed_cm_s = np.sqrt( 8.*8.314*temp*1000./3.14/108. ) * 100.
k_het_N2O5_CLNO2 = Yield_ClNO2 * 0.25 * gamma_N2O5 * ThermalSpeed_cm_s * (AerosolSurfaceArea_um2_cm3*1e-8)
You can find more details in any atmospheric chemistry textbook so I'll cover the bare minimum here. Heterogeneous chemistry is often parameterized as a function of the following parameters/variables:
reactive uptake coefficient (gamma_N2O5): basically the chance of a N2O5 striking a surface that'll trigger a reaction. Dimensionless. 0.03 is a reasonable guestimate here. In reality it depends on a lot of things, such as the nitrate and chloride content on the aerosols, also liquid water content.
Yield_ClNO2: out of 100 N2O5 reactive uptake, 60 ClNO2 is produced... thus 0.6. A reasonable guestimate. In reality this would depend on the chloride content on aerosols as well as liquid water content.
Thermal speed: depends on temperature and the nature of the gas (here 108 is the molecular weight of N2O5).
AerosolSurfaceArea_um2_cm3: exactly what it looks like. Here we treat it as a fixed value for the sake of simplicity. In reality, it can be affected by physical and chemical processes in the atmosphere. 200 um2/cm3 is a reasonable guestimate.
Question: what is the unit of the resulting heterogeneous rate coefficient (k_het_N2O5_CLNO2)?
There are underlying assumptions the way we treated the N2O5-ClNO2 chemistry in this example:
This reaction involves N2O5 and chloride (two reactants), but we write it as a pseudo-1st order reaction, assuming N2O5 is the local limiting reagent. This mostly is not the worst assumption, and also is consistent with most laboratory studies measuring N2O5 reactive uptake coefficient.
Assume gas diffusion is not limiting. See the next section for details.
We here use fixed N2O5 uptake coefficient and ClNO2. In reality these two depend on chemical composition, aerosol mixing state, etc, and can vary a lot! We will not cover this here.
I strongly recommend Stephen Schwartz's work (1986), which is one of my all time favorate! From the mass transport point-of-view, heterogeneous reaction between a gas-phase molecule and a particle consists of multiple steps:
diffusion in the gas-phase, which sets the upper bound of gas-to-particle flux;
interfacial transport: not all molecules colliding the surface will stick to the surface, and the probability of the molecule sticking to the surface over all the molecules striking the surface is often called mass accommodation coefficient, which sets the upper bound of the reactive uptake coefficient.
diffusion in the bulk, which is largely affected by the viscosity of the particle (imagine a particle behaves like water, oily, honey, peanut butter, ... all the way to tar).
chemical reaction within the bulk (including hydration).
Rule of thumb: the slowest step limits the overall rate. Gas diffusion (1) is fast in most cases; but for "very sticky" compounds (with high mass accommodation coefficient), gas diffusion (1) might be limiting. This is often the case for radicals. For very soluble species, interfacial transport (2) may be limiting. For example, glyoxal. For very reactive species, chemical reaction (4) may be faster than the diffusion in the bulk (3), this would lead to quick depletion at the surface and possibly accumulation of product(s) at the surface.
In the case of N2O5 uptake, it reacts quite fast at the surface containing chloride: on average a N2O5 molecule would not penetrate more than a few nm (Gaston et al 2016) into the particle, therefore the overall rate is mostly limited by interfacial transport. However, there are circumstances when N2O5 surface uptake is so fast that it's depleted near the surface, and hence the gas diffusion limits the overall rate. This may be the case if you have high surface area density (e.g. laboratory experiments Lopez-Hilfiker et al 2012), and if yes, you may want to correct for the gas diffusion too.
Don't forget to the time-dependent parameters! You know the trick..
# ========================
# loop over all time steps
# ========================
for i,t in enumerate(tout_s):
# --- set time of day
hour_of_day = (t%86400)/3600.
# --- set time-dependent variables
SZA_deg = (90. + (20.-90.) * diurnal_profile_sine(hour_of_day))
SZA = SZA_deg * np.pi / 180.
MLH_m = 500. + (1500.-500.) * diurnal_profile_sine(hour_of_day)
# --- set time-dependent emissions and depositions
dep_rate_O3_sec = dep_vel_O3_cm_s/(MLH_m*100.)
dep_rate_NO2_sec = dep_vel_NO2_cm_s/(MLH_m*100.)
dep_rate_HNO3_sec = dep_vel_HNO3_cm_s/(MLH_m*100.)
dep_rate_N2O5_sec = dep_vel_N2O5_cm_s/(MLH_m*100.)
This mechanism is pretty small so it takes like 15 seconds to finish a 1-day simulation (output every half an hour) on my end. The following three plots show:
(Left) time series of O3, NO, and NO2. Can you tell me why they look the way they look?
(Middle) time series of N2O5 and ClNO2. This is interesting. At night, N2O5 is an important reservior of reactive nitrogen species, therefore you see the production of N2O5 at night. In the meantime, N2O5 can react with aerosols, producint ClNO2, therefore you see the sharp increase of ClNO2. Also, ClNO2 in this case does not have many nighttime sinks, therefore it'll just accumuate to fairly high levels... until sunrise. ClNO2 will undergo rapid photodissociation, so does N2O5.
(Right) the time series of OH radicals and Cl radicals. OH is the most important oxidant. But Cl, produced from ClNO2 photolysis here, also reaches fairly high concentrations, although it's about 2 orders of magnitude lower than OH. Two reasons that make Cl an important oxidant as well: (i) although it's concentration is usually much lower than OH, it's a lot more reactive than OH radicals! The Cl-oxidation rates of many VOCs are 1-2 orders of magnitude faster than their OH-oxidation rates. (ii) Cl peaks earlier in the day than OH radicals, further increases its importance in the morning.
Not only Cl oxidation is faster than OH oxidation, it leaves unique signatures in the atmosphere too! One thing that can be used to examine if Cl oxidation actually plays a role is to look at certain VOC ratios. Impact of dilution on concentrations would be cancelled out.
i-butane and n-butane have fairly similar OH rate coefficients. Therefore, if OH oxidation dominates, the ratio of i-butane / n-butane wouldn't change much over the course of time. However, their Cl rate coefficients are very different! Therefore if Cl oxidation dominates, the ratio of i-butane / n-butane would change drastically,
i-butane and propane are the other way around: their OH rate coefficients are very different but OH rate coefficients are fairly similar. Therefore if OH oxidation dominates, the ratio of i-butane / propane would vastly change; but if Cl oxidation dominates, this ratio wouldn't change much.
Now here is the fun part: you're the boss of the fake world you just created and you can do whatever you want! For example, to examine the VOC signature with OH oxdiation only, you can simply set the rate coefficients of all CL + VOC reactions to zero. Since we'll only use three VOC species, propane (C3H8), n-butane (NC4H10), and i-butane (IC4H10), we can only turn off these reactions. Remember the utility function that lists all sources and sinks of a given compound?
def sources_sinks(in_spc_name):
print('----------------------------------')
print('%s removal:' % (in_spc_name))
for i,r in enumerate(reactions):
reactants = r.replace(' ','').split('=')[0].split('+')
if in_spc_name.replace(' ','') in reactants:
print('R%d: %s' % (i, r))
print('----------------------------------')
print('%s production:' % (in_spc_name))
for i,r in enumerate(reactions):
products = []
for p in r.replace(' ','').split('=')[1].split('+'):
if '*' in p: products.append(p.split('*')[1]) # in case yield is not 1
else: products.append(p)
if in_spc_name.replace(' ','') in products:
print('R%d: %s' % (i, r))
print('----------------------------------')
sources_sinks('C3H8')
Run this cell and it'll print the following:
----------------------------------
C3H8 removal:
R75: C3H8 + OH = NC3H7O2
R76: C3H8 + OH = IC3H7O2
R267: C3H8 + CL = NC3H7O2
R268: C3H8 + CL = IC3H7O2
----------------------------------
C3H8 production:
----------------------------------
That tells us C3H8 has two OH oxidation pathways, R75 and R76. It also has two CL oxidation pathways, R267 and R268. Do the same for NC4H10 and IC4H10:
----------------------------------
NC4H10 removal:
R103: NC4H10 + OH = NC4H9O2
R104: NC4H10 + OH = SC4H9O2
R269: CL + NC4H10 = NC4H9O2
R270: CL + NC4H10 = SC4H9O2
----------------------------------
NC4H10 production:
----------------------------------
----------------------------------
IC4H10 removal:
R191: IC4H10 + OH = IC4H9O2
R192: IC4H10 + OH = TC4H9O2
R271: CL + IC4H10 = IC4H9O2
R272: CL + IC4H10 = TC4H9O2
----------------------------------
IC4H10 production:
----------------------------------
Now you know what reactions to touch! Let's say if you want to turn off their OH oxidation pathways, in the big loop right before the call to the ODE solver, add these lines:
# --- OH oxidation
rate_coeff[191] = 0.
rate_coeff[192] = 0.
rate_coeff[103] = 0.
rate_coeff[104] = 0.
rate_coeff[75] = 0.
rate_coeff[76] = 0.
Similarly, you can also turn of Cl oxidation and examine only OH oxidation! Do do this, instead of setting those OH rates to zero, you set the Cl rates to zero:
# --- CL oxidation
rate_coeff[271] = 0.
rate_coeff[272] = 0.
rate_coeff[269] = 0.
rate_coeff[270] = 0.
rate_coeff[267] = 0.
rate_coeff[268] = 0.
Toggle back and forth this chunk of codes can turn on/off OH and Cl oxidation pathways! Screenshots below:
Now you need to run the model 3 times
Default, keep both OH oxidation and Cl oxidation on. After the simulation is completed, run this line to archive the concentrations:
conc_both = conc.copy()
Cl oxidation only. Follow the instruction above, turn off OH oxidation. After the simulation is completed, run this line to archive the concentrations:
conc_CL_only = conc.copy()
OH oxidation only. Follow the instruction above, turn off CL oxidation. After the simulation is completed, run this line to archive the concentrations:
conc_OH_only = conc.copy()
After those three sets of simulations, you should have three Pandas DataFrames: conc_both,conc_CL_only,conc_OH_only. Then run the following scripts to make a plot!
marker_list = ['o','s','P']
label_list = ['both', 'Cl only', 'OH only']
for i,c in enumerate([conc_both,conc_CL_only,conc_OH_only]):
plt.scatter(c['IC4H10']/c['C3H8'], c['IC4H10']/c['NC4H10'],
c=tout_s/3600., s=80, edgecolor='w',linewidth=1,
marker=marker_list[i], cmap='inferno',
label=label_list[i])
plt.colorbar().set_label('hour of day')
plt.ylabel('i-butane/n-butane')
plt.xlabel('i-butane/propane')
plt.legend()
Cool huh? As you can see that with OH only, the ratio-ratio line evolves sort of horizontally. But with Cl only, the line evolves vertically! And if you have them both on, the line is kinda line a mixture of both!
In this chapter we briefly introduced heterogeneous, and we examined simple N2O5 heterogeneous reaction on chloride-containg particles producing ClNO2, which could be an important source of Cl radicals. We also examined the unique VOC fingerprint of Cl oxidation in the atmosphere.
One trick we learnt is that, we can run different sensitivity tests by shutting off certain chemical reactions! This is quite useful.
Last but not least, chlorine chemistry in the atmosphere is way more complicated than we had in this example. For example, chlorine chemistry is the major culprit of the ozone hole. Chlorine chemistry also plays a role in the Arctic boundary layer.