[BUG/ISSUE] sea salt alkalinity and heterogeneous acid-catalyzed halogen chemistry

beckyalexander commented 1 year ago

What institution are you from?

University of Washington

Description of the problem

Sea-salt alkalinity was not appropriately considered in the calculation of heterogeneous acid-catalyzed reactions of halogens on sea salt aerosols. The computations were being performed assuming that 100% of sea salt aerosol alkalinity had been titrated. In other words, aerosol pH calculations in ISORROPIA were being applied to total fine-mode sea salt aerosol, and all coarse-mode sea salt aerosol was assumed to have a pH of 5, no matter how much of the alkalinity was being titrated.

Description of troubleshooting performed

This fix will mainly affect simulations where sea salt debromination is turned on. It is currently turned off by default. The attached difference plots show the comparison with sea salt debromination turned on in both the dev and ref simulations.

GEOS-Chem version

13.4.1

Description of modifications

I have fixed this by calculating two new variables in fullchem_HetStateFuncs, f_Alk_SSA and f_Alk_SSC, in KPP that compute the fraction of fine and coarse mode sea salt aerosol that is alkaline. This fraction is then used in fullchem_RateLawFuncs when calculating the first-order uptake of gas onto acidic sea salt aerosol. I multiply (1 - f_Alk_SSA) by the sea salt surface area so that only the acid fraction of the sea salt surface area is used in the rate calculations. I consider this a bug fix because this is how the calculation was performed before the sulfur and sea salt alkalinity calculations were moved into KPP. BrClI_hetchem.txt

Figures

I'm showing difference plots for BrO and O3 for July 2019. I can provide more species upon request. The figure below shows that this fix causes BrO to decrease everywhere as expected. O3 is more complicated, showing increases in the SH and decreases in the NH. ssdebromBrClI_SROSO4s_O3BrO.pdf

I expected increases in O3 everywhere. The reason for the increase in O3 in the NH is entirely due to heterogeneous iodine chemistry. The figure below shows the same difference plots except that the dev simulation only fixed the Br and Cl chemistry (not the I chemistry). So it is the fix to the iodine chemistry that is causing the decrease in NH ozone. The iodine chemistry needs more looking into and may be related to some heterogeneous iodine chemistry that was removed in the more recent versions of GEOS-Chem. I think this is more of a science issue than a bug fix so I'll consider that separately. alkssdebromBrCl_alkssdebrom_BrOO3.pdf

yantosca commented 1 year ago

Thanks @beckyalexander. I'll create a PR and test it. We should be able to get this into 14.1.0, as it is a bug fix.

djxjacob commented 1 year ago

Thanks Becky!

Turning off SSA debromination in the standard model to get ozone right has been a crutch that we would like to get rid of, but the effect from your bug fix seems relatively small. Also the NH decrease is counter to what we want – and your plot without the iodine fix also shows a NH decrease? Clearly the bug fix is important, but it doesn’t seem like it can justify turning SSA debromination back on in the standard model. Let me know if you disagree!

If we keep SSA debromination on as an option, then Bob’s suggestion of bringing your bug fix into 14.1 seems best. If we were to turn SSA debromination back on in the standard model that would require a separate benchmark.

Ccing Randall and Mat in case they are not on this email list.

Daniel

beckyalexander commented 1 year ago

I agree with you @djxjacob that this is not enough to justify turning on sea salt debromination in the standard model, unfortunately. I did notice some Iy sinks that were recently removed from KPP (https://github.com/geoschem/geos-chem/pull/923/files). I put them back in and am running that simulation now. I'm considering that a separate issue though because they are not acid-catalyzed reactions and they were intentionally removed, so it may not be considered a bug.

djxjacob commented 1 year ago

Iodine chemistry is going to change a lot (maybe in 14.2?) with the JPL update to IO + IO. Kelvin has plots on this. Daniel

msulprizio commented 1 year ago

The fix for this (https://github.com/geoschem/geos-chem/pull/1548) has now been merged into the dev/14.1.0 branch.

msulprizio commented 1 year ago

The updates for this fix fail a 1-month benchmark simulation (sea salt debromination off) with KPP convergence errors:

 Forced exit from Rosenbrock due to the following error:
 --> Step size too small: T + 10*H = T or H < Roundoff
 T=   11.964464782736020      and H=   8.4522004874768389E-015
 ### INTEGRATE RETURNED ERROR AT:           24          28           5
===============================================================================
GEOS-Chem ERROR: KPP failed to converge after 2 iterations!
 -> at Do_FullChem (in module GeosCore/FullChem_mod.F90)
===============================================================================

===============================================================================
GEOS-Chem ERROR: Incorrect species units after Do_FullChem!
 -> at Do_Chemistry  (in module GeosCore/chemistry_mod.F90)
===============================================================================

===============================================================================
GEOS-CHEM ERROR: Error encountered in "Do_Chemistry"!
STOP at  -> at GEOS-Chem (in GeosCore/main.F90)
===============================================================================
 Timer forced to stop due to error: GEOS-Chem                     
 Timer forced to stop due to error: All chemistry                 
 Timer forced to stop due to error: => FlexChem

See the full log file here.

@beckyalexander @msl3v @yantosca Any thoughts on how to resolve this as it'a currently holding up 14.1.0 bencharking?

msulprizio commented 1 year ago

@yantosca points out that the log file indicates the issue is likely in the reactions below that return infinity for RCONST:

 ########################################################################
 ### KPP DEBUG OUTPUT!
 ### Reaction rates at problem box           14          13           1
 ########################################################################
                  Infinity SALAAL + O3 + SO2 --> LOx + PSO4 + SO4 - SALAAL
   6.8799391718909510E-009 SALAAL + HCl --> SALACL
   2.3078344094468533E-008 SALAAL + HNO3 --> LOx + NIT
                  Infinity SALCAL + O3 + SO2 --> SO4s + LOx - SALCAL
   5.9665973977075313E-012 SALCAL + HCl --> SALCCL
   9.0869745784508394E-012 SALCAL + HNO3 --> LOx + NITs
   ...

beckyalexander commented 1 year ago

I wonder if a change I made in an earlier patch is in the version that you are running now. In fullchem_SulfurChemFuncs, in the calculation of K_MT(1) and K_MT(4) the below changes were made. Did this get incorporated into this version?

KPP/fullchem/fullchem_SulfurChemFuncs.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/KPP/fullchem/fullchem_SulfurChemFuncs.F90 b/KPP/fullchem/fullchem_SulfurChemFuncs.F90 index 172bf361c..2ec9ea081 100644 --- a/KPP/fullchem/fullchem_SulfurChemFuncs.F90 +++ b/KPP/fullchem/fullchem_SulfurChemFuncs.F90 @@ -204,8 +204,8 @@ CONTAINS srMw = SR_MW(ind_SO2) )

    ! Assume SO2 is limiting, so recompute rxn rate accordingly

K_MT(1) = kIIR1Ltd( C(ind_SO2), C(ind_O3), k_ex ) &
/ C(ind_SALAAL)
K_MT(1) = kIIR1Ltd( C(ind_SO2), C(ind_SALAAL), k_ex ) &
/ C(ind_O3) ENDIF

!------------------------------------------------------------------------ @@ -260,8 +260,8 @@ CONTAINS srMw = SR_MW(ind_SO2) )
```
! Assume SO2 is limiting, so recompute rxn rate accordingly
```
K_MT(4) = kIIR1Ltd( C(ind_SO2), C(ind_O3), k_ex ) &
/ C(ind_SALCAL)
K_MT(4) = kIIR1Ltd( C(ind_SO2), C(ind_SALCAL), k_ex ) &
/ C(ind_O3) ENDIF

!------------------------------------------------------------------------ --

beckyalexander commented 1 year ago

In my view, the + and - signs did not come through here. In each change above, the first two bullet points are minus signs and the second two are plus signs. Since SALAAL and SALCAL are limiting, these should go into the numerator, and O3 should be in the denominator.

yantosca commented 1 year ago

@beckyalexander your changes were added OK.

We are thinking that something in the merge up to KPP 3.0.0 is causing this. I am looking into the issue.

stale[bot] commented 1 year ago

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.

msulprizio commented 1 year ago

The fix for this (#1702) has now been merged into the dev/14.2.0 branch.

msulprizio commented 1 year ago

This update is being removed from 14.2.0 after the initial 1-year full-chemistry GCClassic benchmark simualtion showed drastic decreases in O3 in November and December. See the full discussion here: https://github.com/geoschem/geos-chem/issues/1880. In attempting to debug this, reverting PR #1702 was the only was I was able to resolve the issue.

@beckyalexander has offered to help look into implementing this into a future version of GEOS-Chem, so I will leave some debugging notes here:

I found that the O3 began decreasing first over the Middle East in early November 2019 and spread out from there. I tracked this back to a spike in BrSALA concentrations that began on 11/05/2019 at 22:00 UTC in grid box I=45,J=28, L=9.
This issue also did not show up in a November 2019 run with 4x5 GEOS-FP met. This issue also did not occur in Nadia Columbi's 2016 run with 14.2.0. Both of these suggest that some met condition present in the MERRA-2 fields on this day seems to trigger the bromine explosion over the Red Sea and coastline of Yemen/Saudia Arabia.
Implementing PR #1702 introduced numerous "INTEGRATE RETURNED ERROR" messages from KPP into the GEOS-Chem log file. In the 1-year benchmark simulation, there were ~250 of these error messages in the first spinup month (July 2018) which decreased to ~150 error messages in the final month of the benchmark (Dec 2019). GCClassic didn't fail on these errors which indicates the integrate was successful upon the second attempt. However, the large number of integrate errors shouldn't be ignored -- for reference, we didn't have a single one of these errors in the 14.0.0 1-year benchmark nor in the 1-month benchmark for 14.2.0-alpha.4 (prior to merging in #1702).
The 1-year GCHP benchmark did fail with the integration error on 05/03/2019. We were able to get around this by adding a kludge to adjust the relative tolerance from 0.5e-2 to 0.5e-3 after the first integration failure (see #1878).
In the original updates, RTOL was modified from 0.5e-2 to 0.5e-3 but this caused significant slow down in GEOS-Chem as discussed in #1814. To avoid this slow down, we changed RTOL back to 0.5e-2.

msulprizio commented 1 year ago

This issue should now be addressed by PR #1933.

msulprizio commented 1 year ago

This has now been addressed in #1952

geoschem / geos-chem