NOx chemical loss rate output from fullchem simulation

[cschooling]

Name and Institution (Required)

Name: Chloe Schooling Institution: The University of Edinburgh

Confirm you have reviewed the following documentation

• [x] Support guidelines
• [x] User manuals
• [x] Debugging GEOS-Chem and HEMCO errors
• [x] Current and past Github issues

Description of your issue or question

I am using v 14.0.2.

I would like to output the loss rates of NO and NO2, however this is not an option within the ProdLoss collection. If I want to output this information, would it be best to add some new HISTORY diagnostics, e.g. by modifying Headers/state_diag_mod.F90 (if so, could you direct me on the best way to do this?). Or would this be better/possible to do so by identifying all chemical equations within KPP/fullchem/gckpp_Monitor.F90 where the chemical NO is lost, and determining the loss rate by saving relevant equation rates out to RxnRates collection?

Thank you so much for your help!

[yantosca]

Thanks @cschooling for writing. I think the best way to proceed would be for you to create KPP families for NO and NO2 loss in the fullchem.kpp file as follows:

#FAMILIES                 { Chemical families for prod/loss diagnostic }
POx : O3 + NO2 + 2NO3 + PAN + PPN + MPAN + HNO4 + 3N2O5 + HNO3 + BrO + HOBr + BrNO2 + 2BrNO3 + MPN + ETHLN + MVKN + MCRHN + MCRHNB + PROPNN + R4N2 + PRN1 + PRPN + R4N1 + HONIT + MONITS + MONITU + OLND + OLNN + IHN1 + IHN2 + IHN3 + IHN4 + INPB + INPD + ICN + 2IDN + ITCN + ITHN + ISOPNOO1 + ISOPNOO2 + INO2B + INO2D + INA + IDHNBOO + IDHNDOO1 + IDHNDOO2 + IHPNBOO + IHPNDOO + ICNOO + 2IDNOO + MACRNO2 + ClO + HOCl + ClNO2 + 2ClNO3 + 2Cl2O2 + 2OClO + O + O1D + IO + HOI + IONO + 2IONO2 + 2OIO + 2I2O2 + 3I2O3 + 4I2O4;
LOx : O3 + NO2 + 2NO3 + PAN + PPN + MPAN + HNO4 + 3N2O5 + HNO3 + BrO + HOBr + BrNO2 + 2BrNO3 + MPN + ETHLN + MVKN + MCRHN + MCRHNB + PROPNN + R4N2 + PRN1 + PRPN + R4N1 + HONIT + MONITS + MONITU + OLND + OLNN + IHN1 + IHN2 + IHN3 + IHN4 + INPB + INPD + ICN + 2IDN + ITCN + ITHN + ISOPNOO1 + ISOPNOO2 + INO2B + INO2D + INA + IDHNBOO + IDHNDOO1 + IDHNDOO2 + IHPNBOO + IHPNDOO + ICNOO + 2IDNOO + MACRNO2 + ClO + HOCl + ClNO2 + 2ClNO3 + 2Cl2O2 + 2OClO + O + O1D + IO + HOI + IONO + 2IONO2 + 2OIO + 2I2O2 + 3I2O3 + 4I2O4;
PCO : CO;
LCO : CO;
PSO4 : SO4;
LCH4 : CH4;
PH2O2 : H2O2;
LNO: NO      {add this here, for loss of NO family}
LNO2: NO2    {add this here, for loss of NO2 family}

and then rebuild the KPP solver code:

cd /path/to/your/GCClassic/src/GEOS-Chem/KPP
./build_mechanism.sh fullchem

then recompile and rerun GEOS-Chem. LNO and LNO2 will then be included in the Loss_?LOS? HISTORY ProdLoss collection.

For more information also see:

• The FAMILIES section of kpp.readhedocs.io
• Using KPP with GEOS-Chem

Let us know if you still encounter problems.

NOTE: Normally we would advise you to use the custom mechanism (which is a copy of the fullchem mechanism), but there are some build issues in your version that will be solved in the GEOS-Chem 14.2.1 release (upcoming). So for now it is OK to update the fullchem mechanism. You might also want to keep a backup copy of the fullchem.eqn in case you have to revert to the "out-of-the-box" settings.

[cschooling]

@yantosca Thank you so much for your quick advice.

Creating KPP families makes sense. However for some reason when I try to run ./build_mechanism.sh fullchem/

I get this error: ./build_mechanism.sh: line 99: kpp: command not found KPP failed to build 'gckpp_Rates.F90'! Aborting.

I have tried this also with the original fullchem.kpp file before my changes but it gives the same. Do you know why I might be seeing this error?

Thanks Chloe

[yantosca]

Thanks for writing @cschooling. It looks like you don't have the KPP executable built. You can download and build the KPP executable according to these instructions.

[cschooling]

Thank you

[cschooling]

Sorry didn't mean to close this!

I got that all to work well.

One thing I wondered is how to find the loss rate of NOx in kg/cm3/s. The model gives us the output in molecule/cm3/s, however it is unknown how many of the lost molecules are NO and how many are NO2. Is there a way to find out this information?

Thanks so much!

[yantosca]

Thanks @cschooling. I think you would have to have separate loss families LNOx, LNO, LNO2 to try to back out how much NO and NO2 is lost. I don't think there is a good way to try to back out e.g. LNO and LNO2 just from LNOx.

You could convert the molec/cm3/s of each of those loss families to kg/cm3/s knowing the molecular weights. In the gckpp_Global there is a MW array which has the MW of each species in g.

[yantosca]

@cschooling: MW (g) of NO is MW(ind_NO); MW (g) of NO2 is MW(ind_NO2). These are the same values as are defined in the species_database.yml file.

I also forgot... you would need to add any new prod/loss families (LNO, LNO2) to the species_database.yml file so that there will be defined molecular weights.

[cschooling]

Thanks! I have been using the RxnRates collection to identify in more detail the individual reactions that are contributing to the losses and products of NOx. However, I have come across some discrepancies where there seem to be individual chemical reactions which report a larger rate than the overall product or loss rate, which they are contributing to. For example reaction index 163 (MPN --> PNO2 + PNOx + MO2 + NO2) often gives a larger rate than the rate of both PNO2 and PNOx. Could you explain why this might be?

[yantosca]

Thanks @cschooling. I would defer to @msl3v who knows the KPP families mechanics better than I do.

[msl3v]

Hello @cschooling

Just one quick question: all the P/L species are in units of concentration, and represent absolute change over a timestep. To get the rate from a P/L species, you have to be sure to divide by the timestep. Are you doing this? If not, do that first. If so, then could you provide your gckpp_Monitor.F90 file?

[cschooling]

Hi @msl3v thanks for your help on this!

I haven't divided by the timestep, however the units of the product and loss rate of NO/NO2/NOx (I take from the ProdLoss collection) appear to be in molec cm-3 s-1, which is the same as the units for the reaction rates (I take from the RxnRates collection). I have attached one of my ProdLoss collection .nc4 files, and a screenshot of the headers showing the units. Perhaps I am misunderstanding something?

Thanks so much!

GEOSChem.ProdLoss.20190101_0000z.nc4.zip

[msl3v]

Ah. Thanks @cschooling I hadn't realized that the Prod/Loss for NOx, etc were part of the GEOS-Chem diagnostics. So, I guess we can say that the units are correct.

Let's see... The reaction rates output by the RxnRates diagnostics are the actual reaction rates used by the integrator, e.g. A + B -> C has the rate k[A][B] while the rates output by the Prod/Loss diagnostics are the actual fluxes through the reactions themselves.

Is it possible that at a given timestep & scenario, the rate may be large, but the reactants titrate out? Thus limiting the flux through a reaction at a given rate?

e.g. A = 10, B = 100 (k=1) Rate = 1000, but the max flux can only approach 10?

[cschooling]

Ah ok, makes sense! @msl3v

I am interested in which specific reactions are contributing most significantly to the loss and products of NOx/NO2/NO. Is there a way to determine which specific reactions contributed to the actual fluxes?

[msl3v]

There is!

There are two ways (that I'm aware of) to do this. Both require some code surgery.

Option 1: add a dummy species to each reaction in question. e.g. in fullchem.eqn O3 + NO2 = O2 + NO3 : GCARR_ac(1.20d-13, -2450.0d0); you would add a reactant, say R1, so that O3 + NO2 = O2 + NO3 + R1: GCARR_ac(1.20d-13, -2450.0d0); ... and repeat this for each reaction. Notes:

1. You'll add a different R{X} value for each reaction, e.g. R1... R5.
2. Add the species to the #DEFVAR list at the top of fullchem.eqn (e.g. R1 = IGNORE;)
3. Rerun KPP to generate the necessary *.F90 (if you get to this point, you'll need to have KPP installed. There's good documentation for this.)
4. ... route the species through GEOS-Chem so it can be output as a diagnostic. I'm not sure the procedure here. @yantosca, can you point to directions for this?

If this seems complicated, it's definitely doable.

Option 2: Not so different from (1) You download KPP, check out a branch that has 'flux' enabled. It has a switch that will automatically add dummy species to ALL of the reactions. Then reprocess the fullchem mechanism with this switch 'on'. It will vastly increase the compute burden of running KPP, and the memory required, etc. It really only removes the need to do steps (1) & (2) above by hand. You'll still have to manage (3) & (4). Yet, if you really really really want to do some deep dives on reaction mechanisms & their behavior, this is the way to go. This is super powerful.

[msl3v]

... note: Using reaction fluxes like this, if a reactant has a stochiometric coeff > 1, you'll have to account for that in the reaction flux value.

e.g. A + 2B => C + R1 {R1 being the dummy species} the net change in B is 2R1.

... in case this wasn't obvious..

[cschooling]

@msl3v Thanks so much! I have added dummy species to all the reactions in fullchem.eqn (I call these RA, RB, RC, RD, RE, RF) that I want to monitor. I tried also amending the fullchem.kpp file so I could track their respective loss and products: i.e. adding LRA : RA; (as I did before with NO, NO2, NOx)

I have then rebuilt the KPP solver, and recompiled GEOS-Chem. However, if I try to run the model it is not happy and I get the error attached.

I assume I am missing a step in defining these dummy species so they can be output in diagnostics? @yantosca Do you know how I might do this?

Thanks so much!

[msl3v]

@cschooling, just to be sure, Are you executing the command: $ kpp gckpp.kpp

[yantosca]

@msl3v: The command would be

$ cd GCClassic/src/GEOS-Chem/KPP
$ ./build_mechanism fullchem

which is just a convenience wrapper for kpp gckpp.kpp.

[yantosca]

@cschooling: Make sure also to add the metadata for the new prod/loss species to species_database.yml. It should look similar to that of LOx:

LOx:
  FullName: Dummy species to track loss rate of Ox
  Is_Gas: true
  MW_g: 48.00

It needs to have a molecular weight defined. Also you have to tag it with Is_Gas: True so that GEOS-Chem will know it's a gas species. That is a sanity check -- if a species isn't tagged as either a gas or aerosol, GEOS-Chem will stop the run.

[cschooling]

Thank you @yantosca. Yes I have executed the ./build_mechanism fullchem command. And I have now added the new prod/loss species to the species_database.yml file. However, I still see the same error as above when I run the model..?

[cschooling]

I have realised that I don't have any loss values of these dummy species, so I have removed these loss species. However, weirdly I now get an error: GEOS-Chem ERROR: Could not locate KPP prod/loss species: LCH4! I haven't edited anything related to LCH4 that I'm aware of - so not sure why this has come up.

I have attached the species_database.yml file I use, as well as the fullchem.eqn file and the gckpp_Monitor.F90 file. Archive.zip

[stale[bot]]

This issue has been automatically marked as stale because it has not had recent activity. If there are no updates within 7 days it will be closed. You can add the "never stale" tag to prevent the Stale bot from closing this issue.