Embedded Files
OVERVIEW
OVERVIEW
We now know how to manually set up the chemical mechanism. If you wanna quickly trying things out, it's quite handy. But you might want to avoid this if you want to look at more complicated chemistry!
In this Chapter, we will load the KPP file generated from the MCM website, which is a very powerful tool. Full MCM has like tens of thousands of reactions. But let's start from simple: we will play with CH4 chemistry now (29 species & 68 reactions, including basic inorganic chemistry).
We will also try to turn on the light to get the photochemistry going!
Hit the BROWSE tab, scroll down and find methane in the alkane section. Then click the "+" button. Then scroll all the way up, you should see your mark list, which has only CH4 in it.
Then hit the EXPORT tab, make sure the "include inorganic reactions" and "include generic rate coefficients" are checked.
Then select KPP as the output format, and click the DOWNLOAD button.
You should now have a file called "mcm_export.eqn", which is the KPP file contains all the CH4 chemistry (including all products from CH4 oxidation all the way to CO), and of course basic inorganic reactions (such as NOx, O3, etc). Open that in your favorite text editor and you can see what it looks like.
For the sake of this tutorial, let's rename this file to "mcm_export_ch4.eqn".
p.s. this file is also provided in the GitHub resposiroty (link)
2. Load KPP into the box model
2. Load KPP into the box model
Download this all-in-one Jupyter Notebook (YAMCHA_tutorial_ch2) from GitHub.
Open Google Colab and sign in, then upload the Jupyter Notebook. P.s. You may access GitHub directly if you have a GitHub account and are familiar with GitHub.
Check the screenshots on the right: click that small file folder icon as indicated by the small red arrow. Then drag and drop your KPP file there!
Note: You will get a warning message, reminding you that if you close this runtime, you'll lose this file. That's fine for now. You can choose to mount Colab with your Google Drive (but clicking the button with a small Google Drive icon on it), so you can store your data onto your Google Drive and access them from your Colab. For this tutorial, uploading files into the temporary storage will be just fine.
At this point you should have a file named "mcm_export_ch4.eqn" in your file folder. Right click on it, and select "Copy path". You'll then get the local path (in Colab) which you will be needing later.
We went thru all major cells in Chapter 1 already. Run the first 3 cells, which import libraries, and load a bunch of utility functions that are needed.
Cell #4 is important (which was commented out in Chapter 1 but we will use it here, so make sure it's uncommented). These two lines loads the KPP file. Make sure you replace the file path with the one you got above. In my case it's just "/content/mcm_export_ch4.eqn"
ChemMech = load_mcm_kpp('/content/mcm_export_ch4.eqn')
ChemMech.T
After running that Cell #4 you should get a big table displayed right under that cell, showing the list of reactions and the list of rate coefficients. There's another RO2 list and I'll get back to that later.
And in Chapter 1 we manually edited a chemical mechanism called "ChemMech_simple". We don't need that in this Chapter. So you can comment out the entire Cell #5, or just skip it.
Now go to Cell #6 which is the main workflow of the box model. You're already familiar with it. Let's run the model for 3 days, and output results every 600 seconds (Lines 31 and 31). Also make sure that you set in_ChemMech to ChemMech (Line 38) which is the KPP file you loaded before.
3. Run the model without any sunlight
3. Run the model without any sunlight
Now let's set other user inputs:
# --- this is the solar zenith angle for the j-value parameterization (MCM)
SZA = 90. * np.pi / 180.
# --- environmental conditions
temp = 298 # Kelvin
press = 101325. # Pa
RH = 60.
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['CH4'] = (1800.) *press*(6.0232E+8)/(8.314*temp)
conc_init['NO'] = (10.) *press*(6.0232E+8)/(8.314*temp)
conc_init['O3'] = (30.) *press*(6.0232E+8)/(8.314*temp)
As you can see here, solar zenith angle (SZA) is set to 90 degrees for now, which means there's pretty much no sunlight. Temperature and pressure are 298 K and 101325 Pa. Relative humidity (RH) is set to 60%.
Then the initial concentrations - let's set CH4 to 1800 ppb, NO and O3 to 10 ppb and 30 ppb, respectively.
Then run this cell!
Now take a look at the model outputs. In the very beginning, O3 rapidly drops by 10 ppb, and so does NO. In the meantime, NO2 increased by 10 ppb. This is because of the very rapid NO titration that we mentioned in Chapter 1.
But after that NO2 started to decrease gradually. In fact if you check other nitrogen-containing species (e.g. NO3, N2O5, HNO3), these remain low as well. Since mass has to be conserved - where does all the nitrogen go?
If you plot conc['NA'], i.e. a species called NA, you'll see that it becomes the major reservoir of reactive nitrogen. After 1 day of simulation, it even becomes the dominant. I also have this utility function that can prints all production and removal pathways of a given species. If you run this function, you'll see that NA is produced from R32 and R33 but there's no removal pathway for this thing. So what exactly is this NA?
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('NA')
This is a MCM thing: because HNO3 and N2O5 are the end products of NOx oxidation in the gas-phase, if you keep the model running (without deposition, multiphase chemistry) these species will accumulate to pretty high concentrations. To avoid this, MCM folks set a semi-arbitary lifetime for HNO3 and N2O5 and convert them into "NA". I'm guessing NA stands for "nitrate aerosol".
Similarly, MCM also has a species called "SA".
NOTE: If you run MCM for your own study, you wanna be careful with NA and SA. If you have deposition (we will talk about this later) and / or heterogeneous chemistry or multiphase chemistry (we will talk about this later too), you're probably already getting rid of these gas-phase end products, therefore you should turn off the production of NA and SA.
PRO TIP: How to shut off reactions we don't like?
Find the reaction you want to shut off. In this demo, by running the sources_sinks('NA') you'll find that R32 and R33 are the two production pathways of NA.
In the main loop, right before the call to the ODE solver, insert these two lines (highlighted in yellow) and set their rate coefficients to zero.
# --- now solve the ODE. Use BDF since it's mostly stiff
# accuracy controlled by relative/absolute tolerance (rtol,atol)
# integrate step by step, since some env variables may change over the course of time
rate_coeff[32] = 0. # NA production. This is a MCM thing...
rate_coeff[33] = 0. # NA production. This is a MCM thing...
if i==0: conc[i,:] = conc_init
else: conc[i,:] = solve_ivp(dcodt,tout_s[0:2],conc[i-1,:],t_eval=[output_frequency_seconds],
jac=jac, method='BDF', rtol=1e-5,atol=1e+2).y[:,0]
4. Let There Be Light
4. Let There Be Light
To turn on the light in the model, simply set j-values to non-zero.
We're using photolysis parameterization from MCM, which is a function of solar zenith angle (SZA). Click me for details.
In the model, the big FOR loop, find the line where SZA is, then set it to 30 degrees.
# ==============
# User inputs...
# ==============
# --- this is the solar zenith angle for the j-value parameterization (MCM)
SZA = 30. * np.pi / 180.
The term (Pi/180) converts degree into radius, which is the unit required by the Python built-in trigonometric functions. This is probably vaguely relevant for mid-latitude noontime condition in summer.
Run the model again. You'll see the results look quite different! You actually make quite some O3 and peak OH concentration reaches ~2.5e+7 molec/cm3.
5. Rate/Budget Analysis
5. Rate/Budget Analysis
To answer some important science question, a rate/budget analysis can be quite useful. For instance, in this mechanism (with CH4 chemistry), what are the major OH production/removal pathways? We'll explore this a bit to demonstrate how this can be done.
Firstly, we need to go thru the list of reactions and find all the OH production/removal reactions. You can use this utility function:
sources_sinks('OH')
this should print the following:
----------------------------------
OH removal:
R13: O3 + OH = HO2
R14: H2 + OH = HO2
R15: CO + OH = HO2
R16: H2O2 + OH = HO2
R18: HO2 + OH = PROD
R20: NO + OH = HONO
R21: NO2 + OH = HNO3
R22: NO3 + OH = HO2 + NO2
R25: HO2NO2 + OH = NO2
R27: HONO + OH = NO2
R28: HNO3 + OH = NO3
R30: OH + SO2 = HSO3
R45: CH4 + OH = CH3O2
R56: CH3NO3 + OH = HCHO + NO2
R58: CH3OOH + OH = CH3O2
R59: CH3OOH + OH = HCHO + OH
R63: HCHO + OH = CO + HO2
R67: CH3OH + OH = HCHO + HO2
----------------------------------
OH production:
R12: O1D = OH + OH
R17: HO2 + O3 = OH
R23: HO2 + NO = NO2 + OH
R26: HO2 + NO3 = NO2 + OH
R37: H2O2 = OH + OH
R41: HONO = NO + OH
R42: HNO3 = NO2 + OH
R59: CH3OOH + OH = HCHO + OH
R60: CH3OOH = CH3O + OH
----------------------------------
Then scroll down to the last cell which runs the budget analysis and makes a plot. I'll take OH production as the example and break it down for you. You can find the full cell with codes in the Jupyter Notebook
Obtain the lists for all OH production, then calculate the total OH production rate. Basically add all the sources together.
# --- get the lists of all the OH production/loss
OH_prod_list = ['R12','R17','R23','R26','R37','R41','R42','R59','R60']
# --- calculate the toal OH production rate
total_OH_prod = np.zeros(len(tout_s))
for r in OH_prod_list: total_OH_prod = total_OH_prod + reaction_rates[r]
Then we setup a threshhold of 10% (cutoff_frac = 0.1), and any pathway contributes >10% total OH production or removal, we call it "important". To calculate this we need to integrate the rates over time. We'll use cumulative_trapezoid in scipy.integrate.
# --- select important OH production pathways
cutoff_frac = 0.1 # contributions larger than this is considered as being important
Then we setup the plot, and plot only the "important" pathways, and add those "unimportant" ones to another array and call it "all others".
for r in OH_prod_list:
# --- integrate the rates
integ_total_OH_prod = cumulative_trapezoid(y=total_OH_prod,x=tout_s)[-1]
integ_OH_prod = cumulative_trapezoid(y=reaction_rates[r],x=tout_s)[-1]
if (integ_OH_prod/integ_total_OH_prod)>=cutoff_frac: # plot the pathways with contributions larger than the threshhold
axs[0].plot(tout_s,plotting_scaling_factor*reaction_rates[r],label=r+': '+in_ChemMech.T['reaction'][r])
else: # all the rest go to "others"
others_OH_prod = others_OH_prod + reaction_rates[r]
Finally plot the total, also the "all others" (i.e. the sum of all the unimportant).
You should get the plot on the right. The upper half shows the productions, with the thickest line representing the "total" OH production. The lower half shows the removals, with the thickest line representing the "total" OH removal.
You can see that the most important OH production is HO2 + NO = NO2 + OH, i.e. OH produced from HO2. Because OH and HO2 are rapidly cycling between each other (as a result OH + HO2 is often defined as HOx), this is commonly considered as the secondary OH source. Another pretty important one is the O1D = OH + OH (H2O is also involved here; two OH radicals produced from this reaction). This is the primary OH production, and O1D is produced from O3 photolysis.
As for OH removal, CH4 + OH is a pretty important one, also HCHO + OH. Generally speaking, VOCs are pretty important OH sinks, and in this simple CH4 mechanism we don't have a lot of VOCs for sure. Also NO2 + OH = HNO3 is also pretty important in the beginning, until NO2 is all consumed. Apparently for OH removal, the sum of the all the "unimportant" ones actually plays a pretty big role, so you can reduce the threshold and check more details.
IMPORTANT: Look at the codes producing this plot, there is actually a big bug. Can you spot it?
6. Summary
6. Summary
Now you know how to load KPP files from MCM, and you may generate whatever mechanism you want off MCM website.
Note that if you try to run complicated chemistry (e.g. full MCM), it may take quite some time. You might adjust total run time and output frequency accordingly. The ODE solvers in Python Scipy is pretty okay all-purpose solver, but might not be as fast as some solvers specifically designed for this type of problems.
Nevertheless, in my own work I regularly run mechanisms with nearly 10k reactions, I have also run multiphase chemistry with explicit phase-transfer (naturally these are very stiff & computationally expensive), and I have a version of multiphase 1-dimensional model (with >10k fully coupled ordinary differential equations) using this 0-D model as a building block. This model surely can handle some tough tasks.
Most importantly, more complicated does not necessarily mean more correct. Think before you act. Clearly define your research question(s) then design experiments accordingly. Allow me to share something with you:
"...Models are not right or wrong; they’re always wrong. They’re always approximations. The question you have to ask is whether a model tells you more information than you would have had otherwise. If it does, it’s skillful…"
-- Dr. Gavin A. Schmidt (NASA Goddard)
Google Sites
Report abuse