Interim fix for runaway HMS chemistry

[sdeastham]

Name and Institution (Required)

Name: Sebastian Eastham
Institution: Massachusetts Institute of Technology

Confirm you have reviewed the following documentation

• Contributing guidelines

Describe the update

This update makes two connected changes, both related to the reaction HMS + OH + SO2 -> 2SO4 + CH2O as discussed by Moch et al. (2020) (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020JD032706) which itself implements a collapsed form of the chemistry discussed in Jacob (1986) (https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/JD091iD09p09807). Under the existing implementation, this reaction is implemented as

HMS + OH -> 2SO4 + CH2O - SO2

with a rate constant calculated as k (K_CLD(6) in the code). This assumes that SO$_2$ is always in excess and is not limiting, resulting in SO$_2$ being depleted but its abundance not affecting the rate of reaction. However, in certain circumstances - especially at high spatial resolution - this assumption can become invalid, resulting in integrator failure (see #1982 where @msl3v has performed an analysis and confirmed this problem in standalone testing).

This pull request implements a temporary fix, while a more robust solution is sought. First, the KPP mechanism now represents the reaction as

HMS + OH + SO2 -> 2SO4 + CH2O

and thus removes the issue of a negative product. However, this would result in the reaction being too fast by a factor of [SO2], given that the reaction rate constant k is unchanged. Accordingly, the rate constant is now also changed to be k/[SO2]. This fix is imperfect because the rate constant is set prior to the integrator call and will not update as SO$_2$ concentrations change between integrator time steps. The most likely result is that the overall rate of reaction will now be underestimated in most places, as small reductions in SO$_2$ will result (incorrectly) in the overall rate of reaction slowing down. However, depletion of SO$_2$ during the integrator time step cannot now result in the kind of error previously observed.

Tagging the working group chairs for aerosol (@beckyalexander, @theloniuspunk, @wporter, Fangqun Yu not on GitHub?) and chemistry (@barronh, @luhu0, @jingqiumao) for visibility and @msl3v due to his invaluable diagnostic work and ongoing role in finding a more permanent solution.

Expected changes

This approach has been tested in two settings. First, testing of the bug fix in the setting of #1982 (a C90 simulation) has shown that it succeeds in preventing the integrator failure observed there. Second, three one month tests at C24 have been performed in an otherwise unchanged fullchem simulation. Three situations are tested:

• "old": the original HMS chemistry.
• "new": the new approach.
• "disabled": a simulation in which the reaction is disabled entirely (k set to zero).

A comparison of the monthly mean surface concentration is given below. It shows that the new approach does indeed reduce the rate of loss of HMS and SO2, and subsequently decreases formation of SO4. However the overall change is substantially less than that which results from disabling the reaction entirely.

In terms of the global means surface concentrations, the changes are as follows. This suggests that the overall reduction in HMS processing rate is around 1/3 of what could be expected if the reaction were disabled entirely.

Species	New - Old	Disabled - Old
HMS	+1.9%	+6.8%
SO2	-0.013%	-0.016%
SO4	-0.14%	-0.48%

This picture is further reflected in zonal plots of the monthly mean change.

Reference(s)

Moch et al. (2020) (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020JD032706)
Jacob (1986) (https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/JD091iD09p09807)

Related Github Issue(s)

#1982

[lizziel]

Also tagging @christophkeller since relevant to GEOS-Chem in GEOS.

[lizziel]

@msulprizio, is this going into 14.3.0? @sdeastham, do you need it merged now?

[sdeastham]

@lizziel - I'd say it's moderately urgent? I've had intermittent failures with GCHP at resolutions other than C24, and they've ended up being quite expensive because they are difficult to predict (for example - lost 3 years worth of C180 simulation results when the integrator died during the last month rendering the rest of the results questionable). My only reluctance is that the kludge here, while it does seem to work, is pretty clunky. I'm certain @msl3v would be able to design something cleaner (and/or we could look into something like a limited Update_RCONST along with a two-regime calculation of k_cld(6)) but I'm nervous about the possibility of more lost compute.

[msulprizio]

Yes @lizziel @sdeastham. I will working on bringing this into dev/14.3.0 asap.