Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[BUG/ISSUE] Extremely low mean OH in Tropchem simulation versus Benchmark (v13.0.0 & v13.0.2) #732

Closed
Twize opened this issue May 20, 2021 · 26 comments
Assignees
Labels
category: Bug Something isn't working
Milestone

Comments

@Twize
Copy link

Twize commented May 20, 2021

Description of the problem

I am interested in running troposphere-only full-chem simulations with GEOS-Chem v13. However, for all of the simulations I have run thus far (both on our local machines here at the University of Toronto, and on AWS), I have found the mass-weighted mean OH to be less than half of that from the benchmark (~5 x 10^5 molec. cm-3 versus ~11 x 10^5 molec. cm-3 for the benchmark). The photolysis rates for O(1D) and O(3P) are also extremely reduced in the trop-chem simulation versus the benchmark. This issue leads to a large over-estimation in CO columns, as well as other species including O3.

Description of troubleshooting performed

I have tried multiple trop-chem simulations on both our local machine, as well as on AWS (using the GC v13.0.2 tutorial AMI as a template), and with various settings (i.e., MERRA-2 and GEOS-FP meteorology, 47 and 72 model levels, online and offline emissions) as well as the default un-modified settings when creating a run directory. In all cases, the mass-weighted mean OH given by metrics.py is less than half than that of the benchmark. I have also run benchmark simulations on both the UofT system and on AWS, and comparisons with the official benchmark show very minimal differences, and approximately the same mean OH. I have included a plot comparing the OH fields from the 1-month AWS simulations below:

OH_rst_compare_20190801_Default_vs_benchmark_AWS

I have compared the J-values for O(1D) and O(3P) for the default trop-chem simulation and benchmark I ran on AWS. The photolysis rates for O(1D) are approximately an order of magnitude lower in the tropchem simulation (using default settings), and several orders of magnitude lower for O(3P). See comparison plots below:

O3O1D_Jvals_compare_20190716_AWS_benchmark_vs_tropchem
O3O3P_Jvals_compare_20190716_AWS_benchmark_vs_tropchem

At this point, I am unsure whether it is a bug in the v13 GEOS-Chem codes, or with the the "default" settings that are provided when setting up a run directory, but the differences are apparent on both the systems I have tested, and I have run out of ideas as to what may be the cause. Any help would be appreciated!

GEOS-Chem version

version 13.0.0 (UofT system, git commit: b6f7aff)
version 13.0.2 (AWS, git commit: b36f8efcebc667e932a1d358ed497c1e88c87ded)

The AMI (which is based upon the 13.0.2 tutorial AMI) that I used for the AWS simulations is:
ami-015aec90f03d8f2d3

Description of modifications

No modifications.

Log files

Files below are for the tropchem simulation I ran on AWS (with default settings), had to change the file handle to .txt for some to be able to upload them.

Software versions

  • CMake version:
    3.19.5 (UofT system)
    3.19.2 (AWS)

  • Compilers (Intel or GNU, and version):
    Intel compiler v19.0.3.199 (UofT system)
    GNU Fortran compiler v13.0.2 (AWS)

  • NetCDF version:
    4.2 (UofT system)
    4.7.4 (AWS)

@Twize Twize added the category: Bug Something isn't working label May 20, 2021
@yantosca
Copy link
Contributor

Thanks for writing @Twize. In 13.0.0 there was an issue that met fields were not being read in properly (see #700, #686, #688). I am not sure if one of these is causing your issue.

The photolysis rates for O1D and O3P are I believe computed differently depending on whether or not you are using the full-atmosphere or troposphere only. That might have some impact on your results.

I am tagging @msulprizio, who worked on the implementation of the trop-only UCX mechanism in GEOS-Chem, who may have more insight.

@Twize
Copy link
Author

Twize commented May 25, 2021

Thanks for writing @Twize. In 13.0.0 there was an issue that met fields were not being read in properly (see #700, #686, #688). I am not sure if one of these is causing your issue.

The photolysis rates for O1D and O3P are I believe computed differently depending on whether or not you are using the full-atmosphere or troposphere only. That might have some impact on your results.

I am tagging @msulprizio, who worked on the implementation of the trop-only UCX mechanism in GEOS-Chem, who may have more insight.

Hi Bob, thanks for the quick reply! That bug you mention may possibly be contributing to the issue on our UofT machines, since I am running v13.0.0 there, but it was my understanding that the issue was remedied with v13.0.2 (which is the version I am running on AWS), so don't think it is the cause in that case.

I could always try doing 1-month of a full-atmosphere simulation on AWS to see if the issue with the OH fields persists. I tried this on the UofT machines, and the difference in the mean mass-weighted OH was quite small between the simulations (still approximately half of the benchmark).

Hopefully @msulprizio may have an inkling on what might be the cause of these issues!

@msulprizio
Copy link
Contributor

Hi @Twize. I see you provided the GEOS-Chem dry run log file above. Could you also upload the GEOS-Chem log file for your production simulation(s)? I'm wondering if the restart file is not being found/read correctly.

@Twize
Copy link
Author

Twize commented May 25, 2021

Hi @Twize. I see you provided the GEOS-Chem dry run log file above. Could you also upload the GEOS-Chem log file for your production simulation(s)? I'm wondering if the restart file is not being found/read correctly.

Hi @msulprizio, it seems that I didn't output a log file when I ran the actual simulation. I can re-run the simulation and write the output to a log file (it only takes a couple of hours on AWS). I will upload a log file when its done!

@Twize
Copy link
Author

Twize commented May 26, 2021

Hi again @msulprizio, I am attaching the log file from the 4x5 1-month trop-chem run I did on AWS. Please let me know if there are any other files that you might need, or tests that I can run
that may be useful in diagnosing the issue!

GC.log

@yantosca
Copy link
Contributor

@Twize: It does look like the met data are being read in properly -- you can see that the met fields are being read in to the end of the month (i.e. look for "Found all A1 met fields", etc.)

It also looks like the restart file input is good (as you can see that the values are mostly in the ppt to ppb range for most species (except for NRO2, NAP, SO4{H1,H2,H3,H4}). To me, that indicates the restart file has been read OK.

@Twize
Copy link
Author

Twize commented May 27, 2021

@Twize: It does look like the met data are being read in properly -- you can see that the met fields are being read in to the end of the month (i.e. look for "Found all A1 met fields", etc.)

It also looks like the restart file input is good (as you can see that the values are mostly in the ppt to ppb range for most species (except for NRO2, NAP, SO4{H1,H2,H3,H4}). To me, that indicates the restart file has been read OK.

@yantosca that's good, at least that hopefully rules out the input/met files as the cause of this issue. Let me know if there is anything else that you can think of that I can test which might give us more leads on the issue.

@Twize
Copy link
Author

Twize commented May 28, 2021

@yantosca @msulprizio, I ran another simulation on AWS, this time with troposphere + stratosphere (instead of trop-only), and the mean mass-weighted OH is much closer to the benchmark simulation in this run (~10.7 x 10^5 molec. cm-3 vs. ~4.9 x 10^5 molec. cm-3 in the trop-only simulation). So, it would appear that the issue is coming from something to do with the troposphere-only simulation in v13.

I am attaching relevant files for this new run below:
input.geos.txt
HEMCO.log
GC.log
log.dryrun.txt

I plan to plot up some of the fields from the restart file later today and compare against the benchmark and trop-only run, so I can upload those when they are ready.

@Twize
Copy link
Author

Twize commented May 28, 2021

Hi @yantosca @msulprizio, just following up with some plots of the differences in the restart files, as well as the J-values.

In terms of the OH fields, here is the Benchmark (which I ran myself on AWS) compared with the strat+trop chem run (also on AWS): OH_rst_compare_20190801_AWS_benchmark_vs_fullchem

The OH concentrations are relatively close and where there are differences (max. difference on the order of 0.08 ppt), they are quite inhomogeneous. Contrast this with the Benchmark versus the trop-chem only run (first figure I included in my original post). In that figure, the bias in the OH concentrations is strongly negative almost everywhere, with a maximum difference on the order of ~0.3 ppt.


Plots of the O(1D) and O(3P) photolysis rates tell a similar story, with the benchmark and strat+trop chem simulation matching somewhat closely (at least on the same order of magnitude):

O3O1D_Jvals_compare_20190716_AWS_benchmark_vs_fullchem
O3O3P_Jvals_compare_20190716_AWS_benchmark_vs_fullchem

Contrast these with the second and third figures from the original post, where the photolysis rates were significantly lower in the trop-chem only simulation, particularly for O(3P).


Lastly, to illustrate the impact of the OH differences, I plotted the CO fields from the restart files at the end of the simulations.
Here is the Benchmark versus the strat+trop chem run:

CO_rst_compare_20190801_AWS_tropchem_vs_fullchem_(fixed_scale)
and then the benchmark versus the troposphere-only run:
CO_rst_compare_20190801_AWS_benchmark_vs_tropchem_(fixed_scale)

The positive bias in the CO concentrations is significantly higher compared to the benchmark in the trop-chem only simulation versus the strat+trop chem simulation.

In short, the differences between the benchmark and the strat+trop chem run versus the benchmark and the trop-chem only run are quite drastic, particularly in the photolysis rates and the mean OH. To me, this maybe points to an issue with the trop-chem only simulation is run in GC v13, but maybe you two might have a better idea on what exactly might be the cause of this?

@stale
Copy link

stale bot commented Jun 27, 2021

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.

@stale stale bot added the stale No recent activity on this issue label Jun 27, 2021
@Twize Twize changed the title [BUG/ISSUE] Extremely low mean OH in Tropchem simulation versus Benchmark (v13.0.0 & v13.0.2) [BUG/ISSUE] Extremely low mean OH in Tropchem simulation versus Benchmark (v13.0.0 & v13.0.2) [Never stale] Jun 28, 2021
@stale stale bot removed the stale No recent activity on this issue label Jun 28, 2021
@Twize Twize changed the title [BUG/ISSUE] Extremely low mean OH in Tropchem simulation versus Benchmark (v13.0.0 & v13.0.2) [Never stale] [BUG/ISSUE] Extremely low mean OH in Tropchem simulation versus Benchmark (v13.0.0 & v13.0.2) Jun 28, 2021
@stale
Copy link

stale bot commented Jul 28, 2021

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.

@stale stale bot added the stale No recent activity on this issue label Jul 28, 2021
@r-pound
Copy link
Contributor

r-pound commented Jul 30, 2021

I've been able to replicate this issue in version 13.1.0 with the model built with chemistry only in the troposphere (Tropchem OH average of ~4x10^5). Rebuilding the model with chemistry in troposphere and stratosphere (Standard OH average of ~10x10^5) resolves this so it does appear to be a Tropchem only issue

@stale stale bot removed the stale No recent activity on this issue label Jul 30, 2021
@mje1398
Copy link

mje1398 commented Jul 30, 2021

@yantosca It would seem to me that we should pull TropChem? From what I remember it's no faster than standard and creates a problem as we don't benchmark it? If that's the case I don't really see an advantage to it?

Copy link
Contributor

@mje1398 I believe in the 13.x series, we only have a single mechanism (fullchem) which is UCX-based. You can use a trop-only version of this when you generate a run directory. But as we said, we only benchmark the 72 level version.

@mje1398
Copy link

mje1398 commented Jul 30, 2021

Thanks @yantosca. Perhaps we should remove the option to generate tropchem via the rundirectory as it doesn't work?

@Twize
Copy link
Author

Twize commented Jul 30, 2021

I've been able to replicate this issue in version 13.1.0 with the model built with chemistry only in the troposphere (Tropchem OH average of ~4x10^5). Rebuilding the model with chemistry in troposphere and stratosphere (Standard OH average of ~10x10^5) resolves this so it does appear to be a Tropchem only issue

Glad to know that the issue is reproducible and wasn't only occurring on my end! Thus far, I have switched over to the trop + strat chem version of the model and have not had any issues since in my simulations.

@mje1398
Copy link

mje1398 commented Jul 30, 2021

Glad you've found a solution. I'm thinking that we should just ditch the tropchem version as I don't think its being very useful.

@msulprizio
Copy link
Contributor

It would seem to me that we should pull TropChem? From what I remember it's no faster than standard and creates a problem as we don't benchmark it? If that's the case I don't really see an advantage to it?

About a year ago, the GCST tried to push for retiring the tropchem simulation. The GCSC agreed that we should have a single chemical mechanism but decided they wanted to keep the option for a reduced chemistry grid (see the GCSC meeting minutes with that discussion here). With that said, the tropchem simulation was never properly validated and as indicated by this discussion seems problematic. I suggest we bring up this topic again at the next GCSC meeting next month. I'll contact Randall and Daniel to put it on the agenda.

@mje1398
Copy link

mje1398 commented Jul 30, 2021

Hi @msulprizio,

that’s great. From what I remember there was no time penalty for standard vs tropchem? I which case I can’t see a reason for keeping it?

mat

Copy link
Contributor

@mje1398 Neither can I, but the nested-grid folks always want to use the reduced grid. That's probably why we kept it around.

@msulprizio
Copy link
Contributor

This was the breakdown of timing by component when I ran 1-month benchmarks for the standard (72 and 47 level) and tropchem (47 level) chemical mechanisms in GEOS-Chem 12.8.0:
Screen Shot 2021-07-30 at 12 09 17 PM

@msulprizio
Copy link
Contributor

Hi @Twize @r-pound @mje1398 @djxjacob. I believe I was able to track down the cause of this issue to the following lines in fast_jx_mod.F90:

! Test if the UCX mechanism is being used
IF ( Input_Opt%LUCX ) THEN
!==============================================================
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! %%% FOR MECHANISMS WITH UCX %%%
! %%% (standard, benchmark, *SOA*, marinePOA, aciduptake) %%%
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! SPECIAL TREATMENT FOR H2SO4+hv -> SO2 + 2OH
!
! Only allow photolysis of H2SO4 when gaseous (SDE 04/11/13)
!==============================================================
! Calculate if H2SO4 expected to be gaseous or aqueous
! Only allow photolysis above 6 hPa
! RXN_H2SO4 specifies SO4 + hv -> SO2 + OH + OH
ZPJ(L,RXN_H2SO4,I,J) = ZPJ(L,RXN_H2SO4,I,J) * FRAC
!==============================================================
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! %%% FOR MECHANISMS WITH UCX %%%
! %%% (standard, benchmark, *SOA*, marinePOA, aciduptake) %%%
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! SPECIAL TREATMENT FOR O3+hv -> O+O2 (UCX simulation)
!
! [O1D]ss=J[O3]/(k[H2O]+k[N2]+k[O2])
! SO, THE EFFECTIVE J-VALUE IS J*k[H2O]/(k[H2O]+k[N2]+k[O2])
!
! We don't want to do this if strat-chem is in use, as all
! the intermediate reactions are included - this would be
! double-counting (SDE 04/01/13)
!==============================================================
! Need to subtract O3->O1D from rate
! RXN_O3_1 specifies: O3 + hv -> O2 + O
! RXN_O3_2a specifies: O3 + hv -> O2 + O(1D)
ZPJ(L,RXN_O3_1,I,J) = ZPJ(L,RXN_O3_1,I,J) &
- ZPJ(L,RXN_O3_2a,I,J)
ELSE
!==============================================================
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! %%% FOR MECHANISMS WITHOUT UCX %%%
! %%% (tropchem) %%%
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! Change rate of O(1D)+ N2 to be 3.1e-11 at 298K rather
! than 2.6e-11. The temperature dependence remains the
! same, so the constant changes from 1.8e-11 to 2.14e-11
! according to Heard, pers. comm.,2002. (amf, bmy, 1/7/02)
!==============================================================
! Change the rate of O(1D)+H2O from 2.2e-10 to 1.45e-10*
! exp(89/temp) on the basis of Dunlea and Ravishankara
! 'Measurement of the Rate coefficient for the reaction
! of O(1D) with H2O and re-evaluation of the atmospheric
! OH Production Rate'. One of the RSC Journals
! (mje 4/5/04)
!==============================================================
! Updated from JPL2006, the difference is pretty small.
! (jmao,02/26/2009)
!==============================================================
! Additional update from JPL 10-6 for reaction of
! O3 + hv --> HO2 + OH and O1D + H2:
! Includes calculation of k[O3], where
! k[O3] = J[O3]*1.2e-10/k[O1D], where
! k[O1D] = ([O1D]*[H2]+[O1D]*[H2O])/
! ([O1D]*[H2]+[O1D]*[H2O]+[O1D][N2]+[O1D][O2])
! (bhh, jmao, eam, 7/18/11)
!==============================================================
! Inverse temperature [K-1]
ITEMPK = 1.0_fp / TEMP
! Set species concentrations [molec/m3] ???
C_O2 = 0.2095e+0_fp * NUMDEN
C_N2 = 0.7808e+0_fp * NUMDEN
! Added H2 concentration (bhh, jmao, eam, 7/18/11)
! Seasonal variability of H2 may be important,
! but not included in this update (bhh, jmao, eam, 7/18/11)
C_H2 = 0.5000e-6_fp * NUMDEN
RO1DplH2O = 1.63e-10_fp * EXP( 60.0_fp * ITEMPK ) * C_H2O
RO1DplH2 = 1.2e-10 * C_H2
RO1D = RO1DplH2O &
+ RO1DplH2 &
+ 2.15e-11_fp * EXP( 110.0_fp * ITEMPK ) * C_N2 &
+ 3.30e-11_fp * EXP( 55.0_fp * ITEMPK ) * C_O2
! Prevent div-by-zero
IF ( RO1D > 0.0_fp ) THEN
! RXN_O3_2a specifies: O3 + hv -> O2 + O(1D) #1
ZPJ(L,RXN_O3_2a,I,J) = ZPJ(L,RXN_O3_2a,I,J) * RO1DplH2O / RO1D
! RXN_O3_2b specifies: O3 + hv -> O2 + O(1D) #2
ZPJ(L,RXN_O3_2b,I,J) = ZPJ(L,RXN_O3_2b,I,J) * RO1DplH2 / RO1D
ENDIF
ENDIF

If I comment out the line IF ( Input_Opt%LUCX ) THEN and also comment out the lines between ELSE and ENDIF (i.e. the code previously for non-UCX mechanisms) then I am able to reproduce the OH values of a standard (trop+strat chem) simulation. I am validating this with a 1-month simulation and will report back here with the results.

@djxjacob
Copy link

djxjacob commented Aug 26, 2021 via email

@Twize
Copy link
Author

Twize commented Aug 26, 2021

Thanks for continuing to follow-up on the issue! I am glad that you might have found the bug!

@mje1398
Copy link

mje1398 commented Aug 27, 2021

HI @msulprizio, That's great!

msulprizio added a commit that referenced this issue Sep 2, 2021
…om FAST-JX

Previously in fast_jx_mod.F90 we had RXN_O3_2a reserved for the tropchem
mechanism and reaction O3 + hv --> 2OH and RXN_O3_2b reserved for the
UCX-based standard mechanism and reaction O3 + hv --> O2 + O(1D).The
tropchem mechanism was retired in 13.0.0, and at that time the use of
RXN_O3_2a should have been removed from fast_jx_mod.F90. The continued
handling of this reaction in PHOTRADE_ADJ caused too-low global mean OH in
tropchem simulations and too-high O3 concentrations as originally reported
in #732. This is now resolved
by removing any reference to RXN_O3_2a and any special handling of the
O3 reactions for tropchem simulations from fast_jx_mod.F90. Instead, we
now refer to a single RXN_O3_2 reaction to represent the reaction for
O3 + hv --> O2 + O(1D) now used in all full-chemistry simulations.

Signed-off-by: Melissa Sulprizio <mpayer@seas.harvard.edu>

# Conflicts:
#	GeosCore/hco_interface_gc_mod.F90
@msulprizio
Copy link
Contributor

This bug is now fixed in commit 72f44cd off of GEOS-Chem 13.0.0. This commit has been tagged as 13.0.0+TropchemFix. Users of GEOS-Chem 13.0.0 - 13.2.0 can apply this change to their code using the following steps.

  1. Navigate to the GCClassic/src/GEOS-Chem/ directory. Note: Your path may differ if you cloned the GCClassic wrapper and named it something other than GCClassic.
  2. Do git fetch -p to obtain the latest state of the geoschem/geos-chem repository.
  3. Create a new branch at the current state of your repository for merging the fix into using git checkout -b bugfix/tropchem.
  4. Merge the bug fix into your new branch using git merge 13.0.0+TropchemFix.
  5. Depending on your version, you may be prompted to resolve a conflict in GeosCore/hco_interface_gc_mod.F90. The conflicts should just be spacing differences (see likely conflict code below). I recommend taking the code in the HEAD state of the repository and removing the text , RXN_O3_2a. Don't forget to remove the lines <<<<<<< HEAD and ======= through >>>>>>> 13.0.0+TropchemFix
<<<<<<< HEAD
    USE FAST_JX_MOD,          ONLY : RXN_NO2, RXN_O3_1, RXN_O3_2a
    USE HCO_GeoTools_Mod,     ONLY : HCO_GetSUNCOS
    USE Input_Opt_Mod,        ONLY : OptInput
    USE State_Chm_Mod,        ONLY : ChmState
    USE State_Grid_Mod,       ONLY : GrdState
    USE State_Met_Mod,        ONLY : MetState
=======
    USE FAST_JX_MOD,      ONLY : RXN_NO2, RXN_O3_1
    USE HCO_GeoTools_Mod, ONLY : HCO_GetSUNCOS
    USE Input_Opt_Mod,    ONLY : OptInput
    USE State_Chm_Mod,    ONLY : ChmState
    USE State_Grid_Mod,   ONLY : GrdState
    USE State_Met_Mod,    ONLY : MetState
>>>>>>> 13.0.0+TropchemFix
  1. Commit your merge in either the Git GUI (recommended) or via the command line.
  2. Rebuild GEOS-Chem and re-run your tropchem simulations.

I have tested this fix in 13.2.0 with 1-month benchmark simulations. The model performance statistics and OH metrics are included below.

Screen Shot 2021-09-02 at 9 17 53 PM

1-month Simulations with GEOS-Chem 13 2 0 + Tropchem Bug Fix

This fix does not impact GCHP because GCHP does not support the troposphere-only chemistry option.

NOTE: In GEOS-Chem 13.2.1, we will be retiring the troposphere-only chemistry option following the recommendation of the GEOS-Chem Steering Committee. For more details. please see Github issue #839.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
category: Bug Something isn't working
Projects
None yet
Development

No branches or pull requests

6 participants