Siyuan Wang | Research - Ch.3 Chamber

Ch 3. Environmental Chamber

OVERVIEW

- In previous chapters we learned how to load KPP file (generated from MCM website) into the model. We also briefly looked at the photochemistry of CH4, the simplest volatile organic compound (VOC). There are thousands of chemical species even in the cleanest air, perhaps more. We also learned how to run a simple budget/rate analysis, which can be quite useful especailly for radicals.
- In this chapter, we will take a closer look at one chemical species that is emitted from trees, isoprene (C5H8). The global biogenic emission of isoprene is huge. This thing is very reactive in the air, producing a wide range of compounds. Hundreds of species and thousands of reactions are involved!
- One way to study the atmospheric chemistry is to use environment chambers, because the conditions can be well controlled. In this chapter, we will try to model the OH oxidation of isoprene under high and low NOx conditions, and compare the results to the data reported in Paulot et al Atmos. Chem. Phys. 2009 and Paulot et al Science 2009.

1. Generate the isoprene KPP file from MCM website

Browse or search for isoprene (or C5H8) and mark it. Then EXTRACT the KPP file. Remember to include inorganic reactions, and generic rate coefficients.
Then download the KPP file, and rename it to "mcm_export_isoprene.eqn". I also uploaded this file to GitHub (link).
Download the all-in-one Jupyter Notebook (YAMCHA_tutorial_ch3) from GitHub.
Open Google Colab and sign in, then upload the Jupyter Notebook. You may access GitHub directly if you have a GitHub account and are familiar with GitHub.
Then drag & drop the KPP file into Colab, so that you can access the file in Colab.

2. Set up the model

Run the first 3 cells, importing the libraries and then load all the utility functions.
In Cell #4, load the KPP file (mcm_export_isoprene.eqn) you just got. Remember how to get the file path?

ChemMech = load_mcm_kpp('/content/mcm_export_isoprene.eqn')

ChemMech.T

Run this cell. If the mechanism is loaded successfully, you'll get a big table displayed right under this cell, showing the first 5 rows and the last 5 rows. This mechanism is considerably larger than the CH4 mechanism that we tried in the previous Chapter. The isoprene mechansim has 1944 reactions and 611 species.
Then configure a few more things: let's set the total run time to 0.5 day (12 hours), and output results every 600 seconds. Also make sure the isoprene mechanism you just loaded from the KPP file is used (see the last line of the following code snippet).

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

# Model configuration

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

# --- time stamps

total_run_time_days = 12./24. # how long you want to run the model?

output_frequency_seconds = 600 # how often do you save the output?

tout_s = np.arange(0,total_run_time_days*86400,output_frequency_seconds) # output time stamps, unit: second

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

# Which chemical mechanism to use?

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

in_ChemMech = ChemMech

3. Chamber condition: Paulot et al. ACP (2009)

The chamber experiment (Paulot et al. Atoms. Phys. Chem. 2009) was operated at 296.5 K and with RH < 6%.
Initial concentrations: isoprene: 94 ppbv; NO: 500 ppbv; H2O2: 2.1 ppmv. Here H2O2 was used to generate OH radicals; RH was maintained low to minimize H2O2 loss. NO was used because the major focus of that study was high NOx chemistry, also to suppress O3. Set these initial conditions in the model!

4. Mimic chamber photolysis

Photolysis is a bit tricky. In theory, photolysis frequency (j-value) can be calculated using (spectually resolved) actinic flux, absorption cross section, and quantum yield. But we don't have access to the spectually resolved actinic flux data during that chamber study. If you have this, you may calculate j-values by yourself.
What we can do, in a quick and dirty way, is to use the j-value parameterized by MCM to mimic the photolysis in that chamber. In the paper, the authors reported that the photolysis frequency of H2O2 was 3.1e-6 s^-1. In the model, the photolysis frequency of H2O2 is jval[3], see the jvalue mapping table given in Cell #2. Then in Cell #3 we have:

jval[3] = 1.041E-05*np.cos(SZA)**0.723*np.exp(-1*0.279*np.cos(SZA)**-1)

j-values in MCM is parameterized as a function of SZA (we have tweaked this before). You can try a few SZA values. Apparently when SZA = 64 degrees, jval[3] would be 3.03e-6 s^-1, pretty close to the paper. So we'll go with this photolysis setting in this demo.

Obviously this is not the best way, as the MCM photolysis parameterization is for ambient conditions, while the chamber was using blacklights (276 GE350BL) emitting 300-400 nm, with a maximum at 354 nm.
The full user input section looks like this:

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

# User inputs...

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

# --- this is the solar zenith angle for the j-value parameterization (MCM)

SZA = 64. * np.pi / 180.

# --- environmental conditions

temp = 296.5 # Kelvin

press = 101325. # Pa

RH = 6.

M = press*6.0232E+17/8.314/temp # air density. molec/cm3

N2,O2 = 0.78*M, 0.21*M

H2O = ((RH/100.)*6.1078*np.exp(17.269*(temp-273.3)/temp)*100./8.314/temp)*6.0232e+23/1000000.

RO2 = 0 # total RO2 concentration (molec/cm3). MCM needs this. not really saved

RO2_ts = [] # this is a MCM thing

# --- initial condition. need to convert to molec/cm3

conc_init['C5H8'] = (94.) *press*(6.0232E+8)/(8.314*temp)

conc_init['NO'] = (500.) *press*(6.0232E+8)/(8.314*temp)

conc_init['H2O2'] = (2100.) *press*(6.0232E+8)/(8.314*temp)

Run the model!!!

5. Modeled vs measured:

Moment of truth: now let's compare a few more modeled species to the measured!
In the GitHub repository, I uploaded a CSV file with some observations reported in Paulot et al. (ACP 2009), link. Download that file, drag & drop into Colab.
Then in Cell #7, load the CSV file into a Pandas DataFrame:

paulot2009 = pd.read_csv('/content/Paulot_et_al_ACP_2009.csv')

paulot2009

If loaded successfully, a table will be printed immediately below.

Then in Cell #8, I wrote some script to make a quick plot, evaluating the modeled time series of isoprene, O3, HONO, glycolaldehyde, hydroxyacetone, and HO2NO2 using observations from that CSV file:

nrows, ncols = 3, 2

paulot_spc_list = ['isoprene_Paulot2009','O3_Paulot2009','HONO_Paulot2009','GLYC_Paulot2009_ppb','HACET_Paulot2009_ppb','PNA_Paulot2009']

this_work_list = ['C5H8','O3','HONO','HOCH2CHO','ACETOL','HO2NO2']

display_name = ['isoprene',r'O$_3$','HONO','Glycolaldehyde','hydroxyacetone',r'HO$_2$NO$_2$']

fig, axs = plt.subplots(nrows=nrows,ncols=ncols,figsize=(6,6))

for i in range(nrows):

for j in range(ncols):

n = i*ncols+j

axs[i,j].plot(paulot2009['Time_Paulot2009_min'],paulot2009[paulot_spc_list[n]],'ko',label='Paulot et al 2009')

axs[i,j].plot(tout_s/60.,conc[this_work_list[n]]/(press*(6.0232E+8)/(8.314*temp)),'C3-',label='this model')

axs[i,j].set_ylabel(display_name[n]+' (ppb)')

if i==nrows-1: axs[i,j].set_xlabel('time elapsed (min)')

axs[0,0].legend()

plt.subplots_adjust(top=0.95, bottom=0.15, left=0.15, right=0.95, hspace=0.2, wspace=0.3)

Note that you need to convert the units accordingly. As we discussed before, conc is in molec/cm3, but the data in the CSV is in ppb. In the scripts above, all concentrations (molec/cm3) are converted to mixing ratio (ppb).

You should get a plot that looks similar to the figure on the right!
Usually measurements would have uncertainties. Also our model would have uncertainties too! One important source of model uncertainty is from the kinetics, e.g. uncertainties of measured rate coefficient may span vastly from a few % to hundreds of %. Yield / branching ratio estimated based on measurements would have uncertainties too, which are essentially translated / propagated from individual measurements.
With everything considered, I'd say the agreement is not bad at all!

6. Homework:

In the previous sections, we set up the model for Paulot et al Atmos. Chem. Phys. 2009, which was under high NOx condition. With the model, we can also test what's going on under low NOx condition, and we can compare to Paulot et al Science 2009.
Condition of the Paulot et al Science 2009 paper is kinda similar: initial isoprene 94 ppbv, initial H2O2 1.660 ppmv, and let'ts set initial NO to 0. Tune SZA until you get reasonable decay of isoprene. Then plot the time series of isoprene, IEPOX, and ISOPOOH. I'm attaching my plot on the right for your reference.

7. Further notes:

In this chapter, we set up the model in a way that is only remotely comparable to chamber condition, therefore the comparison to chamber measurements, although not bad, is not apple-to-apple. It's natural to use 0-D model to model chambers, but there are things you should pay close attention to. I'm not an expert in this area so I'll only list a few in my experience:

- Wall loss: all molecules in the chamber are constantly colliding with the walls, and once they collide there is a chance they may stick there. Some molecules are more "sticky" than others, e.g., HNO3, H2O2, NH3. Not only gases, particles undergo wall loss as well! Depending on sizes and properties, the wall loss characteristic time could vary a lot. For example, NH4NO3 wall loss is pretty fast. Usually humidity would speed up the wall loss; and wall loss is generally a bigger concern for smaller chambers. One way to account for the wall loss in the model, is to add this loss mechanism using the characteristic time, which may be estimated based on measurements. Let's say you injected some H2O2 into a particle-free chamber at very low RH, and you watch the decay under dark. Using the data you can probably estimate a "wall loss characteristic time scale". Then you can add this into the model, by adding a first-order loss reaction. Let's again take H2O2 wall loss as an example here.

H2O2 =

rate_coeff = 1/characteristic_time_sec

where characteristic_time_sec is just H2O2's wall loss characteristic time (in second). The reaction is just " H2O2 = ", without any products, becuase we don't care what the product is.

Note that sometimes people would report "wall loss corrected data", or normalize the measurements of a reactive constituent to a non-reactive one. In these cases, you probably don't need to worry too much about wall loss. Bottom line is, you need to know what exactly you're modeling and what exactly you're comparing to.

- Photolysis: in this demo we used the MCM built-in parameterization for photolysis, as a function solar zenith angle. This is for ambient conditions (clear sky), and obviously doesn't work for most chamber conditions. Photolysis is a chamber-specific thing. You need to figure out a proper way to estimate the photolysis in your chamber. Sometimes people would estimate a "background" 1st-order photolytical decay rate for certain species and take it as the effective photolysis frequency of this species for the chamber experiments (watch out for wall loss, OH-propagation, etc). You can load this into the model, and scale the j-values of other species to this species.
- Condensation / evaporation particles: we will not discuss this in detail, but there are equations out there.
- Heterogeneous reactions: we will briefly touch this topic in Chapter 5. One thing unique about chamber heterogeneous reactions is that, you may need to correct for the diffusion limitations, and this is especially the case for big particles and very sticky compounds.

Google Sites

Report abuse