Incorrect sign of exponents in several mercury reaction rates (HgY + Z + M -> HgYZ + M, where Y=Br, OH, Cl and Z=NO2, HO2, BrO, ClO)

[arifein]

Name and Institution (Required)

Name: Ari Feinberg Institution: IQF-CSIC and MIT

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 GEOS-Chem v14.3. I think that the reaction rate formulations for several reactions in the second-step of Hg(0) oxidation have the incorrect sign of their exponents. Specifically, the reactions are (from KPP/Hg/Hg.eqn):

-HgBr    + NO2    = HgBrNO2 :  GCJPLPR_abab(4.3d-30, 5.9d0, 1.2d-10, 1.90d0, 0.6d0);
-HgBr    + HO2    = HgBrHO2 :  GCJPLPR_abab(4.3d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgBr    + ClO    = HgBrClO :  GCJPLPR_abab(4.3d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgBr    + BrO    = HgBrBrO :  GCJPLPR_abab(4.3d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgBr    + OH     = HgBrOH :   GCJPLPR_abab(4.3d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgCl    + NO2    = HgClNO2 :  GCJPLPR_abab(4.3d-30, 5.9d0, 1.2d-10, 1.90d0, 0.6d0);
-HgCl    + HO2    = HgClHO2 :  GCJPLPR_abab(4.3d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgCl    + ClO    = HgClClO :  GCJPLPR_abab(4.3d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgCl    + BrO    = HgClBrO :  GCJPLPR_abab(4.3d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgOH    + NO2    = HgOHNO2 :  GCJPLPR_abab(4.1d-30, 5.9d0, 1.2E-10, 1.90d0, 0.6d0);
-HgOH    + HO2    = HgOHHO2 :  GCJPLPR_abab(4.1d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgOH    + ClO    = HgOHClO :  GCJPLPR_abab(4.1d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgOH    + BrO    = HgOHBrO :  GCJPLPR_abab(4.1d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgOH    + Br     = HgBrOH :   GCJPLPR_abab(4.1d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);
-HgOH    + OH     = HgOHOH :   GCJPLPR_abab(4.1d-30, 5.9d0, 6.9d-11, 2.40d0, 0.6d0);

Table 1a in Shah et al. (2021) shows that the rates should have negative exponents. The listed exponent signs in the table agree with the original references Jiao and Dibble (2017), Wu et al. (2020): 

If we look at the rate law for GCJPLPR_abab in Hg_RateLawFuncs.F90, no change in sign is applied to the exponents:

  FUNCTION GCJPLPR_abab( a1, b1, a2, b2, fv ) RESULT( k )
    !
    ! Third body effect for pressure dependence of rate coefficients.
    ! a1, b1 are the Arrhenius parameters for the lower-limit rate.
    ! a2, b2 are the Arrhenius parameters for the upper-limit rate.
    ! fv     is the falloff curve paramter, (see ATKINSON ET. AL (1992)
    !        J. Phys. Chem. Ref. Data 21, P. 1145). Usually fv = 0.6.
    !
    ! For these reactions, these Arrhenius law terms evaluate to 1:
    !    EXP(c1/T)
    !    EXP(c2/T)
    ! because c1 = c2 = 0.  Therefore we can skip computing these
    ! terms.  Also, fct1 = fct2 = 0, so we will skip computing these
    ! terms as well.  This is more computationally efficient.
    ! (bmy, 06 Jan 2022)
    !
    REAL(dp), INTENT(IN) :: a1,   b1,    a2,    b2,   fv
    REAL(dp)             :: rlow, rhigh, xyrat, blog, fexp, k
    !
    rlow  = a1 * ( K300_OVER_TEMP**b1 ) * NUMDEN
    rhigh = a2 * ( K300_OVER_TEMP**b2 )
    xyrat = rlow / rhigh
    blog  = LOG10( xyrat )
    fexp  = 1.0_dp / ( 1.0_dp + ( blog * blog ) )
    k     = rlow * ( fv**fexp ) / ( 1.0_dp + xyrat )
  END FUNCTION GCJPLPR_abab

Proposed Fix

The signs of the exponents have to be set to negative in KPP/Hg/Hg.eqn:

+HgBr    + NO2    = HgBrNO2 :  GCJPLPR_abab(4.3d-30, -5.9d0, 1.2d-10, -1.90d0, 0.6d0);
+HgBr    + HO2    = HgBrHO2 :  GCJPLPR_abab(4.3d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgBr    + ClO    = HgBrClO :  GCJPLPR_abab(4.3d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgBr    + BrO    = HgBrBrO :  GCJPLPR_abab(4.3d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgBr    + OH     = HgBrOH :   GCJPLPR_abab(4.3d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgCl    + NO2    = HgClNO2 :  GCJPLPR_abab(4.3d-30, -5.9d0, 1.2d-10, -1.90d0, 0.6d0);
+HgCl    + HO2    = HgClHO2 :  GCJPLPR_abab(4.3d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgCl    + ClO    = HgClClO :  GCJPLPR_abab(4.3d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgCl    + BrO    = HgClBrO :  GCJPLPR_abab(4.3d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgOH    + NO2    = HgOHNO2 :  GCJPLPR_abab(4.1d-30, -5.9d0, 1.2E-10, -1.90d0, 0.6d0);
+HgOH    + HO2    = HgOHHO2 :  GCJPLPR_abab(4.1d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgOH    + ClO    = HgOHClO :  GCJPLPR_abab(4.1d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);
+HgOH    + BrO    = HgOHBrO :  GCJPLPR_abab(4.1d-30, -5.9d0, 6.9d-11, -2.40d0, 0.6d0);

Impacts

The incorrect signs will yield discrepancies of a factor of 2 or more for temperatures below 273 K. However, this issue has a negligible impact on the global Hg cycle, as these second-stage reactions are minor. I ran a year-long simulation for 2015 to compare the standard v14.3 with this correction, and there are insignificant changes to TGM, wet deposition, etc (see the benchmark). Nevertheless, if one is interested in specific minor species of Hg(II), this error will affect their formation (e.g., 80% less HgBrNO2 in the polar surface atmosphere, see attached plot).

So in summary, I don’t think this affects any completed simulations results, but it would be good to change in case reactions with similar rate formulations are implemented in the future or if users are interested in specific Hg(II) species.

Full benchmark plots: benchmark_5900_5901.pdf

Tagging @viral211 and mercury WG co-chairs @jennyfisher @yanxuz @hmhorow

[viral211]

Thanks for looking into this @arifein. But I am a little confused. From what I see, the subroutine GCJPLPR_abab( a1, b1, a2, b2, fv ) uses the positive exponents b1, b2 for 300/TEMP. For example, rhigh = a2 * ( K300_OVER_TEMP**b2 ). Wouldn't it be consistent with the negative exponents for TEMP/298 in Table 1?

[arifein]

Hi @viral211, maybe I didn't explain this clearly. Table 1 is correct in your paper, but the current GC code is incorrect in including positive exponents (the first code block that I sent). I think the correct version would be with negative values for b1 and b2, right?

[viral211]

Hi @arifein. What I was trying to say is that GCJPLPR_abab(4.3d-30, 5.9d0, 1.2d-10, 1.90d0, 0.6d0) will calculate $k\infty$ as $1.2 \times 10^{-10} (300/T)^{1.9}$, which is equivalent (approx) to $k\infty=1.2 \times 10^{-10} (T/298)^{-1.9}$ from Table 1 snippet in your original post. So the exponents passed to GCJPLPR_abab need to have their signs reversed (to positive) in this case. Sorry if I am reading this wrong.

[arifein]

Hi @viral211, you're totally right, I mis-read the order of the K300_OVER_TEMP variable. Sorry about that misunderstanding, the reactions are listed correctly in v14.3. Time for the weekend :)

[viral211]

No worries! I have also misread this exact piece of code before. 🤪