Embedded Files
OVERVIEW
OVERVIEW
We now know how to manually set up the chemical mechanism. If you wanna quickly trying things out, this is quite handy. But you might want to avoid this if you want to look at more complicated chemistry! But the idea is quite simple: if you already have the species and reaction lists prepared, you can copy-paste them into the BIG TABLE and define the rate coefficients (kr wave) accordingly. You can also manually edit an automatically loaded chemical mechanism.
In this Chapter, we will load the KPP file generated from the MCM website, which is a very powerful tool. Full MCM has probably 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 hit the "Add selection to the mark list". Then scroll all the way up, you should see your mark list, which has only CH4 in it.
Then hit the EXTRACT tab, check three items: KPP, experimental KPP format, Include inorganic reactions, and Include generic rate coefficients. See screenshot:
Then hit the extract button. This triggers the download and you'll get your own chemical mechanism file in KPP format, with CH4 chemistry (including all products from CH4 oxidation all the way to CO), and of course basic inorganic reactions (such as NOx, O3, etc)
2. Load KPP into IGOR
2. Load KPP into IGOR
Firstly open the blank template (pxp). Then click the button indicated by the red arrow. Then in the pop-out dialog, select the KPP file you just downloaded.
Then click the ODE_ScratchBoard button, to see if there is an ODE_ScratchBoard.ipf pops out. If yes, kill it.
Then click the button (indicated by the red arrow), remember this little thingy? This would generate the ODE_ScratchBoard, based on the KPP file you just loaded.
CTRL+M, open the main procedure window, and set the initial concentrations for CH4, O3 and NO, which are 1800 ppb, 30 ppb, and 10 ppb, respectively.. I have the codes for you. Also see the screenshot.
Concentration_Matrix[0][13] = 1800 *Pressure*(6.0232E+8)/(8.314*TEMP) // CH4 Concentration_Matrix[0][4] = 30 *Pressure*(6.0232E+8)/(8.314*TEMP) // O3 Concentration_Matrix[0][9] = 10 *Pressure*(6.0232E+8)/(8.314*TEMP) // NO
3. Set kr wave
3. Set kr wave
Hit the button indicated by the red arrow. This would bring you a text wave with all rate coefficients processed.
Copy the entire wave. CTRL+M, paste the copied kr wave, also inside the FOR loop, something like this.
(I know this is also kinda frustrating and I'll improve this in the future)
4. Ready! Set! GO!
4. Ready! Set! GO!
Set total run time to 3 days, and time interval to 600 sec, i.e. we ask the model to run for 3 days and output results every 10 min. Then hit the biggest button and run the model!
This mechanism is still very simple so it takes a few seconds to finish. After it's done, select species to preview in the list.
Anything interesting? CH4 remains flat, its oxidation products such as HCHO, CO are all zero. NO (initially 10 ppb) immediately reacts with O3 (initially 30 ppb) producing NO2, then NO2 gradually converted into NO3, N2O5, and eventually NA. Now what the heck is NA???
Run the rate analysis tool: type in the command line:
and you'll see these printed in the history window:
--------------------------------------------------------------- NA removal: --------------------------------------------------------------- NA production: 32: HNO3 = NA 33: N2O5 = NA + NA 33: N2O5 = NA + NA ---------------------------------------------------------------As seen there is no removal pathways for NA, and HNO3, N2O5 eventually all converted into 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-random lifetime for HNO3 and N2O5 and convert them into "NA". I'm guessing NA stands for "nitrate aerosol".
Same idea for SO3, which is irreversibly converted into "SA" (maybe it's "sulfate aerosol").
NOTE: If you run MCM for your own study, you wanna be careful about NA and SA. If you have deposition (we will talk about this later) and / 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 case, the formation of NA and SA are reaction No.32, 33, and 34.
CTRL+M, in the big FOR loop, set kr[32], kr[33], and kr[34] to zero! Alternatively, you can comment these out.
5. LET THERE BE LIGHT
5. 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 this snippet:
// =========================== // Photolysis reactions // j-values parameterizations from MCM // =========================== SZA = 90 * Pi / 180 Jvalue[1] = 6.073E-05*cos(SZA)^1.743*EXP(-1*0.474*sec(SZA)) Jvalue[2] = 4.775E-04*cos(SZA)^0.298*EXP(-1*0.08*sec(SZA)) Jvalue[3] = 1.041E-05*cos(SZA)^0.723*EXP(-1*0.279*sec(SZA)) .....
you see we first set SZA to 90 degrees. The term (Pi/180) converts degree into radius, which is the unit required by the IGOR built-in trigonometric functions. SZA = 90 degree means it's dark. So if you want some light, simply set SZA to <90. Let's do this:
SZA = 20 * Pi / 80
This is probably vaguely relevant for mid-latitude noontime condition.
Now let's run the model again! VIOLA! CH4 is still pretty much flat, but it's only because its lifetime is so long. But if you plot the modeled OH concentration:
Display concentration_matrix[][22] vs time_output // OH is species No.22
The modeled OH radical reached 1e+6 molec/cm3! Hey this looks quite reasonable!
And if you look at HCHO, CO, they both increased! Apparently OH produced with the presence of light oxidized some CH4, produced HCHO and CO.
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 a while to simply generate the ODE_Scratchboard.ipf. It will take even longer to run. The IGOR built-in solver is based on CVODE which is one of the best in the business. But IGOR itself is not designed for heavy duty calculations (we have FORTRAN for that).
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. IGOR 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