[FEATURE REQUEST] Mercury chemistry (reduction and oxidation) diagnostics through ProdLoss for v13.4 #1284

arifein · 2022-06-16T17:06:58Z

With the new updates to the Hg simulation and Hg chemistry being added to KPP, we lost the diagnostics for net and gross Hg0 oxidation from the MercuryChem HISTORY collection (ProdHg2from*). In an earlier version of @viral211 's new chemistry mechanism, there were ProdLoss tracers that could be tracked via KPP and outputted in the ProdLoss HISTORY collection. For example, in his Hg.eqn file:

{OXIDATION BY BR}
Hg0     + Br     = HGBR + RM0pBtM1 :        GCARR2(1.46E-32,1.86,0.0);
HGBR             = Hg0 + RMBtM0 :           GCARR2(1.6E-9,1.86,-7801.0);
HGBR    + Br     = Hg0 + RMBpBtM0 :         3.9E-11;
HGBR    + NO2    = HG0 + RMBpNO2tM0 :       3.0E-12;
HGBR    + NO2    = HGBRNO2 + RMBpNO2tM2 :   GCJPLPR(4.3E-30,5.9,0.0,1.2E-10,1.90,0.0,0.6,0.0,0.0);
HGBR    + HO2    = HGBRHO2 + RMBpHO2tM2 :   GCJPLPR(4.3E-30,5.9,0.0,6.9E-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + ClO    = HGBRCLO + RMBpCLOtM2 :   GCJPLPR(4.3E-30,5.9,0.0,6.9e-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + BrO    = HGBRBRO + RMBpBROtM2 :   GCJPLPR(4.3E-30,5.9,0.0,6.9E-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + OH     = HGBROH + RMBpOHtM2 :     GCJPLPR(4.3E-30,5.9,0.0,6.9E-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + Br     = HGBR2 + RMBpBRtM2 :      3.0E-11;
HGBRO   + CH4    = HGBROH + RMBOpCH4tM2 :   GCARR(4.1E-12,0.0,-856.0);
HGBRO   + CO     = HGBR + RMBOpCOtMB :      GCARR(6.0E-11,0.0,-550.0);
HGBR    + O3     = HGBRO + RMBpO3tM2 :      3.0E-11;

These RM* tracers were then outputted by ProdLoss. However, this didn't seem to make it into the release version. I was wondering if this functionality could be added in a similar way to a newer GEOS-Chem version, or if there was a reason this wasn't added in? This would make it easier to benchmark simulations and evaluate how Hg chemistry changes between versions.

The text was updated successfully, but these errors were encountered:

yantosca · 2022-06-21T14:39:36Z

Thanks @arifein for writing. I think this might have gotten removed at some point by me during testing. I can restore these. Sorry for the bother.

yantosca · 2022-06-21T15:57:02Z

A couple more questions for @arifein @viral211 @msl3v:

Going back to the code that I received from @msl3v, the only prod/loss species that were being used were:

{OXIDATION BY BR}
Hg0     + Br     = HGBR + RM0pBtM1 :        GCARR2(1.46E-32,1.86,0.0);
HGBR             = Hg0 + RMBtM0 :           GCARR2(1.6E-9,1.86,-7801.0);
HGBR    + Br     = Hg0 + RMBpBtM0 :         3.9E-11;
HGBR    + NO2    = HG0 + RMBpNO2tM0 :       3.0E-12;
HGBR    + NO2    = HGBRNO2 + RMBpNO2tM2 :   GCJPLPR(4.3E-30,5.9,0.0,1.2E-10,1.90,0.0,0.6,0.0,0.0);
HGBR    + HO2    = HGBRHO2 + RMBpHO2tM2 :   GCJPLPR(4.3E-30,5.9,0.0,6.9E-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + ClO    = HGBRCLO + RMBpCLOtM2 :   GCJPLPR(4.3E-30,5.9,0.0,6.9e-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + BrO    = HGBRBRO + RMBpBROtM2 :   GCJPLPR(4.3E-30,5.9,0.0,6.9E-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + OH     = HGBROH + RMBpOHtM2 :     GCJPLPR(4.3E-30,5.9,0.0,6.9E-11,2.40,0.0,0.6,0.0,0.0);
HGBR    + Br     = HGBR2 + RMBpBRtM2 :      3.0E-11;
HGBRO   + CH4    = HGBROH + RMBOpCH4tM2 :   GCARR(4.1E-12,0.0,-856.0);
HGBRO   + CO     = HGBR + RMBOpCOtMB :      GCARR(6.0E-11,0.0,-550.0);
HGBR    + O3     = HGBRO + RMBpO3tM2 :      3.0E-11;

... etc ...

{OXIDATION BY OH}
HG0     + OH      = HGOH + RM0pOtM1 :       GCARR2(3.34E-33, 0.0, 43.0);
HGOH              = HG0 + RMOtM0 :          GCARR2(1.22E-9, 0.0, -5720.0);
HGOH    + NO2     = HGOHNO2 + RMOpNO2tM2 :  GCJPLPR(4.1E-30,5.9,0.0,1.2E-10,1.90,0.0,0.6,0.0,0.0);               

... etc ...

HGOH    + O3      = HGOHO + RMOpO3tM2 :     3.0E-11;
HGOHO   + CH4     = HGOHOH + RMOpCH4tM2 :   GCARR(4.1E-12,0.0,-856.0);
HGOHO   + CO      = HGOH + RMOOpCOtMO :     GCARR(6.0E-11,0.0,-550.0);

Are there others that should be attached to e.g. the Cl equations? Could someone help me add these back in?

Also -- do we have metadata for the species database? I think we'd need to add for each of these:

PRM0pBtM1:
  Fullname: name of species
  Is_Prod: true

Thanks!!

yantosca · 2022-06-21T16:08:19Z

Also, if you look at the gckpp_Monitor.F90 output, you see that the reactions are:

  CHARACTER(LEN=100), PARAMETER, DIMENSION(30) :: EQN_NAMES_0 = (/ &
     '   Hg0 + Br --> RM0pBtM1 + PRM0pBtM1 + HgBr                                                         ', & ! index 1
     '       HgBr --> RMBtM0 + PRMBtM0 + Hg0                                                              ', & ! index 2
     '  HgBr + Br --> RMBpBtM0 + PRMBpBtM0 + Hg0                                                          ', & ! index 3
     ' HgBr + NO2 --> RMBpNO2tM0 + PRMBpNO2tM0 + Hg0                                                      ', & ! index 4
     ' HgBr + NO2 --> RMBpNO2tM2 + PRMBpNO2tM2 + HgBrNO2                                                  ', & ! index 5
     ' HgBr + HO2 --> RMBpHO2tM2 + PRMBpHO2tM2 + HgBrHO2                                                  ', & ! index 6
     ' HgBr + ClO --> RMBpCLOtM2 + PRMBpCLOtM2 + HgBrClO                                                  ', & ! index 7
     ' HgBr + BrO --> RMBpBROtM2 + PRMBpBROtM2 + HgBrBrO                                                  ', & ! index 8
     '  HgBr + OH --> RMBpOHtM2 + PRMBpOHtM2 + HgBrOH                                                     ', & ! index 9
     '  HgBr + Br --> RMBpBRtM2 + PRMBpBRtM2 + HgBr2                                                      ', & ! index 10

Is this what would be expected, to have the RM* and the PRM* species both listed as products?

viral211 · 2022-06-21T22:19:02Z

I could add the metadata for these diagnostic species to species_database.yml. However, I am not sure if we need to include all of these P/L diagnostics in the standard code. We could simplify this by having diagnostics for gross oxidation of Hg(0) to Hg(II) from each pathway (Br, OH, and Cl) and a diagnostic for gross reduction of Hg(II) to Hg(0). What do you think @arifein?
@yantosca. The diagnostics for Cl were not added because it is a minor pathway, but they should be for completeness. I can add that too. Also, gckpp_Monitor.F90 has two species because I couldn't get the P/L rates unless I assigned the species into P/L families (the PRM* refer to family names). Maybe there is a way to do this without having separate family names. If so, please let me know.

arifein · 2022-06-21T22:37:59Z

@viral211 I agree with you that most users would not need all of the diagnostics and that the gross oxidation by pathway and gross reduction would be enough!

yantosca · 2022-06-23T22:24:07Z

@arifein: Would you be willing to take the point on this?

arifein · 2022-06-24T14:58:22Z

Hi @yantosca, yes I can look into adding this. I will try thinking about how to simplify the diagnostics as well. @viral211 do you think we would still need to include all of the RM* in the Hg.eqn but then just calculate and output the summary diagnostics?

viral211 · 2022-06-24T16:01:16Z

@arifein we could just define production and loss families in KPP/Hg/gckpp.kpp as follows:

#FAMILIES
LHg0 : Hg0 + HgBr + HgOH + HgCl; {gross Hg2 production}
PHg0 : Hg0 + HgBr + HgOH + HgCl; {gross Hg2 loss}
PHg2Br : HgBrNO2 + HgBrHO2 + HgBrOH + HgBrBrO + HgBrClO + HgBrO + HgBr2; {Hg2 production from Br}
PHg2OH : ... ; {Hg2 production from OH}
PHg2Cl : ... ; {Hg2 production from Cl}

Then run kpp to regenerate the gckpp* files. kpp will add the P/L species to the reactions, which you can see in gckpp_Monitor.F90. This should be sufficient to get the oxidation and reduction diagnostics. We don't need to add all of the RM* species.
Let me know if you run into issues. Thanks for doing this.

yantosca · 2022-06-24T16:07:37Z

Thanks @viral211 and @arifein. I agree, it's better to use the #FAMILIES option of KPP.

stale · 2022-07-25T06:47:56Z

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.

arifein · 2022-08-11T18:20:31Z

Hi all, I have been working on adding the mercury chemistry diagnostics to ProdLoss to validate the v14 Hg simulation. Following @viral211's advice, I added to Hg.kpp:

#FAMILIES                 { Chemical families for prod/loss diagnostic }
PHg2 : HgBrNO2 + HgBrHO2 + HgBrClO + HgBrBrO + HgBr2 + HgBrOH + HgBrO + HgOHNO2 + HgOHHO2 + HgOHClO + HgOHBrO + HgOHOH + HgOHO + HgClNO2 + HgClHO2 + HgClClO + HgClBrO + HgClBr + HgClOH + HgClO + HgCl2 + Hg2ClP + Hg2ORGP + HgBr + HgOH + HgCl; {gross GOM production}
PHg0 : Hg0; {gross Hg0 production}
PHg2Br : HgBrNO2 + HgBrHO2 + HgBrClO + HgBrBrO + HgBr2 + HgBrOH + HgBrO + HgBr; {GOM production from Br}
PHg2OH : HgOHNO2 + HgOHHO2 + HgOHClO + HgOHBrO + HgOHOH + HgOHO + HgOH; {GOM production from OH}
PHg2Cl : HgClNO2 + HgClHO2 + HgClClO + HgClBrO + HgClBr + HgClOH + HgClO + HgCl2 + HgCl; {GOM production from Cl}

So here PHg2 tracks the gross oxidation to HgI + HgII, similar to the previous v12 variable from MercuryChem, Gross_Hg_Ox. PHg0 tracks the gross reduction to Hg0, meaning that the net oxidation (ProdHg2fromHg0) can be calculated as PHg2 - PHg0. I get reasonable global values for these variables after a 1-year simulation: PHg2 = 19149 Mg/yr, PHg0 = 12452 Mg/yr.

I also included the diagnostics for the different 1st-step oxidants (Br, OH, Cl). Note that the inclusion of HgBrOH within the Br family is a bit ambiguous because it can also form through the OH pathway, but from Shah et al. (2021) I understand that the more probable pathway for this compound is through Br. There is always the option for users to create more specific reaction diagnostics, where needed, but I think this should be sufficient for most applications. What do you think?

arifein · 2022-08-11T18:41:47Z

As a technical note, some other code changes where necessary to get these diagnostics working:

After running KPP, I had to add the following missing section to gckpp_Util.F90:

 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! Get_OHreactivity - returns the OH reactivity
! The OH reactivity is defined as the inverse of its lifetime.
! This routine was auto-generated using script OHreact_parser.py.
! Generated on 2021-05-05
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Get_OHreactivity ( CC, RR, OHreact )

! CC - Concentrations of species (local)
  REAL(kind=dp) :: CC(NSPEC)
! RR - reaction rates (local)
  REAL(kind=dp) :: RR(NREACT)
! OHreact - OH reactivity [s-1]
  REAL(kind=dp) :: OHreact

  OHreact = RR(9)*CC(95) + RR(22)*CC(88) + RR(26)*CC(96) + RR(33)*CC(98)

END SUBROUTINE Get_OHreactivity
! End of Get_OHreactivity subroutine
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Because of other updates since v12, around line 1232 in mercury_mod.F90 the loops need to be modified:

       ! Chemical loss of species or families [molec/cm3/s]
       IF ( State_Diag%Archive_Loss ) THEN
          DO S = 1, State_Diag%Map_Loss%nSlots
             KppId = State_Diag%Map_Loss%slot2Id(S)
             State_Diag%Loss(I,J,L,S) = C(KppID) / DT
          ENDDO
       ENDIF

       ! Chemical production of species or families [molec/cm3/s]
       IF ( State_Diag%Archive_Prod ) THEN
          DO S = 1, State_Diag%Map_Prod%nSlots
             KppID = State_Diag%Map_Prod%slot2Id(S)
             State_Diag%Prod(I,J,L,S) = C(KppID) / DT
          ENDDO
       ENDIF

I got an error that SPC_LHg* were not found in the restart file, so I modified the flags in HEMCO_Config.rc from EFYO to EY:
* SPC_ ./GEOSChem.Restart.$YYYY$MM$DD_$HH$MNz.nc4 SpeciesRst_?ALL? $YYYY/$MM/$DD/$HH EY xyz 1 * - 1 1
The outputted units of the ProdLoss variables were incorrectly listed as kg/s, when they are actually molec/cm3/s (see code in 2.). A similar issue was mentioned in [BUG/ISSUE] ProdLoss diagnostic units variable is set to Kg s-1 in NetCDF output #404 for the tropchem simulation, because there was an if statement for IsFullChem in Headers/state_diag_mod.F90 setting all other simulations' ProdLoss units to kg/s. I don't know if this would be correct, but I just removed the if statement. I assume that all ProdLoss variables are calculated in molec/cm3/s:

    ELSE IF ( TRIM( Name_AllCaps ) == 'LOSS' ) THEN
       IF ( IsDesc    ) Desc  = 'Chemical loss of'
       IF ( isRank    ) Rank  = 3
       IF ( isTagged  ) TagId = 'LOS'

       ! NOTE: Units are different depending on simulation, due to historical
       ! baggage.  Maybe clean this up at a later point to use the same units
       ! regardless of simulation type. (bmy, 12/4/17)
       IF ( isUnits   ) THEN
          Units = 'molec cm-3 s-1'
       ENDIF

    ELSE IF ( TRIM( Name_AllCaps ) == 'PROD' ) THEN
       IF ( isDesc    ) Desc  = 'Chemical production of'
       IF ( isRank    ) Rank  = 3
       IF ( isTagged  ) TagId = 'PRD'

       ! NOTE: Units are different depending on simulation, due to historical
       ! baggage.  Maybe clean this up at a later point to use the same units
       ! regardless of simulation type. (bmy, 12/4/17)
       IF ( isUnits   ) THEN
          Units = 'molec cm-3 s-1'
       ENDIF

viral211 · 2022-08-12T16:38:39Z

Thank you @arifein. This looks good to me. I thought you would have had to add the prod and loss species to the species_database.yml file. I guess that's not needed anymore.

arifein · 2022-08-15T14:54:05Z

Thanks Viral! I didn't add them so far to species_database.yml and no errors came up, so I guess this was changed from past versions

yantosca · 2022-08-17T13:26:08Z

Thanks @arifein! We can pull this into the code.

Also, the OHreactivity is taken care of in the 14.0.0 dev branch. I've updated the build_mechanism.sh and OHreact_parser.py scripts for the new KPP 2.5.0+. Very soon we'll have the mechanisms rebuilt with KPP 3.0.0 (once we are sure the journal reviewers on Haipeng's paper don't ask for any KPP revisions...)

These updates address the issues raised in geoschem/geos-chem #1284 Headers/diaglist_mod.F90 - Added IsHg and IsCarbonCycle module variables. (The IsCarbonCycle is future-proofing us for 14.1.0, when that mechanism will be added.) Headers/state_diag_mod.F90 - PROD and LOSS for simulations with KPP are now assigned units of molec cm-3 s-1. PROD and LOSS for other specialty simulations are assigned units of kg s-1 (until we bring these into KPP). GeosCore/mercury_mod.F90 - Updated loops that store in to State_Diag%Prod, following suggestion by Ari Feinberg KPP/Hg/Hg.kpp - Added prod families PHg0, PHg2, PHg2Br, PHg2OH, PHg2Cl to Hg.kpp run/shared/species_database_hg.yml - Added entries for PHg0, PHg2Br, PHg2OH, PHg2_Br, PHg2Cl KPP/Hg/gckpp* - Rebuilt with KPP 2.5.0 Signed-off-by: Bob Yantosca <[email protected]>

yantosca · 2022-08-23T20:57:29Z

We can close this issue as PR #1356 has now been merged into the 14.0.1 development branch.

arifein added the category: Feature Request New feature or request label Jun 16, 2022

yantosca self-assigned this Jun 21, 2022

yantosca added specialty simulations labels Jun 21, 2022

yantosca added category: Bug Something isn't working and removed category: Feature Request New feature or request labels Jun 21, 2022

yantosca added this to the 14.0.1 milestone Jun 21, 2022

arifein mentioned this issue Jun 21, 2022

[BUG/ISSUE] Hg simulation in version 13.4 missing land and ocean emissions #1285

Closed

stale bot added the stale No recent activity on this issue label Jul 25, 2022

yantosca added never stale Never label this issue as stale and removed stale No recent activity on this issue labels Jul 25, 2022

yantosca mentioned this issue Aug 22, 2022

Bug fixes for Prod/Loss diagnostics for Hg simulation -- Closes #1284, #1285 #1356

Merged

yantosca closed this as completed Aug 23, 2022

msulprizio modified the milestones: 14.0.1, 14.0.0 Sep 19, 2022

msulprizio added the topic: Hg or POPs simulations Related to the GEOS-Chem Mercury Simulation label Oct 1, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[FEATURE REQUEST] Mercury chemistry (reduction and oxidation) diagnostics through ProdLoss for v13.4 #1284

[FEATURE REQUEST] Mercury chemistry (reduction and oxidation) diagnostics through ProdLoss for v13.4 #1284

arifein commented Jun 16, 2022

yantosca commented Jun 21, 2022

yantosca commented Jun 21, 2022

yantosca commented Jun 21, 2022

viral211 commented Jun 21, 2022

arifein commented Jun 21, 2022

yantosca commented Jun 23, 2022

arifein commented Jun 24, 2022

viral211 commented Jun 24, 2022

yantosca commented Jun 24, 2022

stale bot commented Jul 25, 2022

arifein commented Aug 11, 2022

arifein commented Aug 11, 2022

viral211 commented Aug 12, 2022

arifein commented Aug 15, 2022

yantosca commented Aug 17, 2022

yantosca commented Aug 23, 2022

[FEATURE REQUEST] Mercury chemistry (reduction and oxidation) diagnostics through ProdLoss for v13.4 #1284

[FEATURE REQUEST] Mercury chemistry (reduction and oxidation) diagnostics through ProdLoss for v13.4 #1284

Comments

arifein commented Jun 16, 2022

yantosca commented Jun 21, 2022

yantosca commented Jun 21, 2022

yantosca commented Jun 21, 2022

viral211 commented Jun 21, 2022

arifein commented Jun 21, 2022

yantosca commented Jun 23, 2022

arifein commented Jun 24, 2022

viral211 commented Jun 24, 2022

yantosca commented Jun 24, 2022

stale bot commented Jul 25, 2022

arifein commented Aug 11, 2022

arifein commented Aug 11, 2022

viral211 commented Aug 12, 2022

arifein commented Aug 15, 2022

yantosca commented Aug 17, 2022

yantosca commented Aug 23, 2022