[BUG/ISSUE] Unexpected differences in full-chemistry simulation caused by TOMAS updates

[msulprizio]

Name and Institution (Required)

Name: Melissa Sulprizio
Institution: Harvard/GCST

Confirm you have reviewed the following documentation

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

Description of your issue or question

Following the merge of the PRs listed below into dev/14.2.0, we observed differences in GCClassic and GCHP 1-hour full-chemistry benchmark simulations when we expected none.

• Implementation of TOMAS in GCHP #1378
• Bug fix for TOMAS sulfate production rates #1569
• Implementation of TOMAS in GCHP HEMCO#164
• Implementation of TOMAS in GCHP GCHP#252

Differences are summarized below:

• GCC_Difference_Summary.txt
• GCHP_Difference_Summary.txt

We will revert dev/14.2.0 to before the merge of the TOMAS PRs because we need to better understand these differences.

Tagging @BettyCroft @yantosca

[yantosca]

Upon further investigation, the differences w/r/t the benchmark simulation are due to the following modifications in the fullchem.eqn file:

PH2SO4     = IGNORE; {SO4 from gas-phase chemistry}
PSO4AQ     = IGNORE; {SO4 from cloud chemistry}

//
// Cloud
// S(IV) --> S(VI)
SO2 + H2O2          = SO4 + PSO4AQ       : K_CLD(1);
SO2 + O3            = SO4 + PSO4AQ       : K_CLD(2) + SRO3; {Jan 2023; Added SRO3 here; BA}
SO2 {+O2}           = SO4 + PSO4AQ       : K_CLD(3);        {Mn & Fe catalysis + HET_DROP_CHEM()}

HMS + OH            = 2SO4 + CH2O - SO2 + 2PSO4AQ : K_CLD(6);        {Sep 2021; Moch2020; MSL}

SO2 + OH {+M} = SO4 + HO2 + PH2SO4 :         GCJPLPR_aba(3.30d-31, 4.3d+00, 1.6d-12, 0.6d0);

CH2OO + SO2 = CH2O + SO4 + PH2SO4:           3.70d-11;                                                                         {2019/11/06; Bates2019; KHB}

CH3CHOO + SO2 = ALD2 + SO4 + PH2SO4 :        7.00d-14;                                                                         {2015/09/25; Millet2015; DBM,EAM}

Question: is there a reason the extra species weren't defined under the #FAMILIES section of fullchem.kpp? KPP does some special handling to make sure that family species aren't counted in reactions as active products.

Tagging @BettyCroft @Jourdan-He @msulprizio @lizziel

[yantosca]

@BettyCroft: I think I remember this... it doesn't look like there is a clean way of making the PH2SO4 and PSO4AQ species work as #FAMILIES, as I think KPP doesn't have a way of knowing that in those rxns SO4 is aqueous. I think I remember that is why we proceeded this way.

If we need to keep this setup, then we will have to implement a separate mechanisms for TOMAS (i.e. in its own GEOS-Chem/KPP/tomas subfolder). Then we would have to compile that with -DMECH=tomas -DTOMAS=y -DTOMAS_BINS=15 or 30. Because this causes changes to the fullchem simulation (and the benchmarks), we cannot bring it in as it is now. Happy to consider alternatives.

[yantosca]

Also tagging @theloniuspunk

[BettyCroft]

@yantosca - yes I am remembering this in the same as you are in your previous message. I'm not aware of any alternatives.

@msl3v helped me with these code development to pass the aqueous sulfate production rates from KPP to TOMAS - and perhaps he can weigh in here regarding if there are any alternatives we should consider. My understanding from our previous discussions was that these dummy species were the best option for pasing the cloud chem production rates to TOMAS..

[BettyCroft]

@theloniuspunk, how does the approach of a separate chemical mechanism for TOMAS sound to you?

[theloniuspunk]

I guess this is fine. I worry a bit about how this may affect the TOMAS simulations down the road when the full chem mechanism gets updated. Will these updates now not transfer to the TOMAS simulations?

[BettyCroft]

@yantosca I think that @theloniuspunk brings up an important question here regarding a separate chemical mechanism for TOMAS. Disconnecting TOMAS from future chemical mechanism updates and benchmarking seems like an undesirable option - please send along some clarification about how that would proceed in the future if there was a separate mechanism for TOMAS.

In case this is of help here - I'll note that I was getting floating point exceptions and negative concentrations from a couple of sections of code not directly related to TOMAS in 1) fullchem_HetSateFuncs.F90 and 2) mixing_mod.F90 and I added a couple of lines of code to prevent those issues - could changes coming from either of those subroutines be causing the differences?

[yantosca]

Thanks @BettyCroft. We think that we can use this mechanism as-is but it should be defined as a separate mechanism from the fullchem. We have to be sensitive to not only the numerical differences that are generated w/r/t fullchem. but also adding extra species adds extra memory overhead that can affect how GEOS-Chem works in e.g. the NASA/GEOS model, where memory usage is already at a premium.

NOTE: We are about to freeze 14.2.0 so this might have to be one of the first updates in 14.3.0.

[msl3v]

I'm unable to comment on this in github. Here's my comment: Hey all. Adding any new species to KPP will alter the results, as we've learned. Esp with dummy species like were added here, since they are initialized to zero, their impact on the error norm is not small. This will affect step size, and integration outcomes. With regard to the species themselves, they were added manually, and not as FAMILY species. There's several ways this could be resolved. I'm not sure where to look for the actual results to see the magnitudes of the differences.

[BettyCroft]

@yantosca, would you be able to point out where/how to view the magnitude of the differences?

[msl3v]

You can set the dummy species value to a large number before integration, and then subtract the number afterwards and effectively eliminate the impact on integration. It's awkward but would be a good test.

In a branch of KPP, there is code that will skip over FAMILY P/L species in the error norm calc. It would be possible to enable manually adding species to that list.

But honestly, neither of the options sounds particularly elegant. If we can determine that the non-zero diff is reasonable, I'd just say go with it.

[yantosca]

Interesting idea @msl3v.

[yantosca]

@BettyCroft: The differences showed up in our 1hr benchmark output:

• GCHP vs GCHP OH metrics
• GCC vs GCC OH metrics

[BettyCroft]

Thanks @yantosca and @msl3v, - are these differences low enough to be in the realm of just going with it?

[msl3v]

It would be worth performing the test to make sure this is the cause.

…

On 4/11/23 18:39, BettyCroft wrote: Thanks @yantosca <https://github.com/yantosca> and @msl3v <https://github.com/msl3v>, - are these differences low enough to be in the realm of just going with it? — Reply to this email directly, view it on GitHub <#1742 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVE236UM55XH7Y7XJOK2BLXAXMSBANCNFSM6AAAAAAWWYSQ5Q>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[stale]

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.

[BettyCroft]

Hi @yantosca - checking in on this issue - what are your suggestions on how to proceed? Please let me know if you have suggestions on how I might be able to help. Thanks, Betty

[msulprizio]

Hi @BettyCroft. The plan is to include the TOMAS updates linked above first thing in 14.3.0 and accept the numerical differences as noise from the solver. However, issue #1741, which documents the failure of TOMAS simulations in both GCClassic and GCHP, still remains. We are hoping the TOMAS simulation can be fixed in a separate pull request for inclusion in 14.3.0 or soon thereafter, we just don't want the code to diverge too much from the original pull requests.

[BettyCroft]

Thanks @msulprizio! That would be great if the TOMAS update for KPP (linked above) can go in 14.3.0. In regard to the other updates for implementation of TOMAS in GCHP and removal of TOMAS' dependency on bpch diagnostics - I agree that we do not want the code to diverge too much from the original pull requests. Please let me know if I can be of help.

[yantosca]

I think we can close this issue now that PR #2060 has been merged.