Ch 1(a). Nighttime NOx
OVERVIEW
- In this chapter, let's play around with simple nighttime NOx chemistry!
- Chemical mechanism is the heart of a chemistry model. It's nothing but two lists: (i) all chemical species involved in this system, and (ii) all chemical reactions that may occur. You will also feed the model with these two lists, along with rate coefficients and initial concentrations.
- To provide the model with a chemical mechanism: (i) you can generate & download KPP file from the MCM website; or (ii) you can type-in these by hand! We will cover (i) in the next chapter. In this chapter we will do a demo with the 2nd way.
- Be careful when you type by hand! If you made a typo then you may mess it up! Sometimes IGOR may complain, but sometimes the model would run just fine and mistakes can be extremely difficult to spot.
- This tutorial seems long (because I blah a lot) but it actually takes less than 20 minutes.
- NOTE: in this example, all rate coefficients are constant. If this is the case, everything can be done in the GUI and you do not have to touch the source codes!! So the model can be "coding-free" ;)
STEP-BY-STEP
1. Download the blank template, an IGOR pxp file. Open it! (if for some reason this link is broken, please email me and I'll send you the file)
2. You'll see this GUI. Depending on the version, this GUI might look different. If you don't see it, hit WINDOWS - OtherWindows - iChamber (yes I know it's lame give me a break!)
3. Hit the button BIG TABLE, this brings up a table with 4 (empty) waves*. kr_value is for the rate coefficients, in (molec/cm3)^(1-m) s^-1, where m is the reaction order. All rate coefficients passed into the solver must be in this unit. Other waves are exactly what they look like.
* for folks new to IGOR, the IGOR-folks refer to arrays as waves. The waves can be numerical or other types, and can be 1-dimensional (column), 2-dimensional (row*column), 3-dimensional (row*column*chunk), ...
4. Our simple nighttime NOx chemistry consists of five species (NO, NO2, NO3, N2O5, and O3). Please type these species and reactions into the table. In the Species wave, there should not be spaces. Don't use weird characters in species names. Also type in the initial concentrations in ppb, but in the model everything is converted into molec/cm3 (298K, 1atm).
Note: your species are numbered! In this case, NO is species No.0, NO2 is species No.1, ... and O3 is species No.4. So are your reactions.
5. Enter the rate coefficients. These are the 298 K values.
6. Now let's enter the reactions. A few rules:*
- Add a space before and after each species;
- ... unless the stoichiometric yield is not one, in which case the yield and the species is separated by a * and no spaces in between.
See the example on the right. Each arrow is a space.
* NOTE: The newest version of the model (updated on 1 Dec 2019) will automatically fix the formatting! So you don't have to type of spaces if you aren't a big fan of these rules. Still caution should be taken when manually entering chemical mechanisms.
7. Please hit the button indicated by a red arrow (alternatively click WINDOWS - ProcedureWindows), to see if there is an items called "ODE_ScratchBoard.ipf" pops out. If yes, kill it!
~ UNDER THE HOOD ~
ODE_ScratchBoard.ipf is the procedure contains all the ordinary differential equations. Someone else (a.k.a. me) made this one, but you need to make your own! So let's kill the already existing one.
NOTE: whenever your chemical mechanism is changed, i.e. you added / removed / edited a species / reaction, you need to kill a pre-existing ODE_ScratchBoard and make a new one. If you only changed initial concentration or rate coefficient, you don't have to do this.
8. Now hit this button! This makes your own ODE_ScratchBoard, based on the chemical mechanism you provided. If you have a big mechanism (e.g. hundreds / thousands of species / reactions), this may take a while. Don't disturb! Go get a coffee or pet your cat.
Once it's completed, you can click ODE_ScratchBoard button to view. Don't kill it now (hide the window if you don't wanna see it)! Otherwise you'll have to re-make it.
9. Now again please click WINDOWS - Procedure Windows - ODE_ScratchBoard.ipf, to open your own:
- YDOT is the array that carries the ordinary differential equation for each species, i.e. dC/dt.
- YDOT[n]: n is the species index. Remember that NO is species No.0. So YDOT[0] is really d[NO]/dt. The unit of YDOT[n] is molec/cm3/s.
- kr[i]: rate coefficient of reaction i. In this example, reaction No.0 is NO + O3 = NO2, and reaction No.1 is NO2 + O3 = NO3, etc.
- Conc[k] is the concentration for species k, in molec/cm3.
Let's take a break! Take a look at your reaction list: NO is involved in two reactions:
R0 NO + O3 = NO2
R2 NO + NO3 = 2*NO2
In both R0 and R2, NO is a reactant. So if you write down the rate expression of R0:
- d[NO]/dt = - d[O3]/dt = d[NO2]/dt = k0*[NO]*[O3]
- d[NO]/dt = - d[NO3]/dt = 2*d[NO2]/dt = k2*[NO]*[NO3]
Therefore, total d[NO]/dt = -k0*[NO]*[O3] - k2*[NO]*[NO3]. Total d[NO]/dt here is YDOT[0], k0 and k2 are just kr[0] and kr[2], respectively. [NO] is Conc[0] because NO is species No.0, [O3] is Conc[4] because O3 is species No.4, etc. As you can see, this translates the ordinary differential equations into the format that IGOR can understand. This is all you need to know at this point.
NOTE: this Conc[] wave contains local values at each integration step! This is useless to you! I'll explain more later.
10. We're almost there! Let's go back to the GUI, in the red box, set Total run time (day) to 0.02, and time interval (sec) to 5. This means, we run the model for 0.02 day (about half an hour), and ask the model to output the results every 5 seconds.
11. Then hit the big button with a black triangle to run the model! Watch the line evolves!
This mechanism is fairly simple, so it takes only a few seconds. Once it's finished, click on different species to see the time series.
(Assume you're somewhat familiar with IGOR)
The time wave in the model is called time_output (unit: second). If you type this in the command line:
edit time_output
you'll see this wave. We ran the model for 0.02 day, and output results every 5 seconds, so the number of points of this wave should be:
0.02 * 86400 / 5
The concentrations of each species at each output time step (not integrate step) is stored in a 2-dimensional wave called concentration_matrix (unit: molec/cm3). If you type this in the command line:
edit concentration_matrix
The rows represent time step, i.e. it has the same number of rows as the time_output. The columns represent each species, i.e. the number of columns is equal to the number of species. Therefore, if you wanna plot the time series of NO (which is species No.0), you should type in the command line:
display concentration_matrix[][0] vs time_output
And if you wanna append NO2 and O3 on it, type this:
appendtograph concentration_matrix[][1] vs time_output appendtograph concentration_matrix[][4] vs time_output
Now you should get the left plot, which is ugly! I made the right plot which is slightly less ugly.
So this is what how the system evolves: initially you have some NO, some NO2, and some O3. NO rapidly reacts with O3, producing NO2. This is often referred to as "NO titration". The amount of NO2 increase is quite close to (but not equal to) NO and O3 decrease. If you plot NO3 and N2O5, you'll see that a small amount of NOx (NO+NO2) is converted into NO3 and N2O5. So it all makes sense!
Practice: could you check what total NOy (NO+NO2+NO3+2*N2O5) looks like?
This is it! You now have a 0D box model with simple nighttime NOx chemistry!