How to fix CH4 lifetime in GEOS Chem CH4 simulation？

XunZhangNJU commented 1 month ago

Your name

XunZhang

Your affiliation

Nanjing University

Please provide a clear and concise description of your question or discussion topic.

In the GEOS-Chem CH4 simulation, within the carbon_gases_mod.F90 file, the following code exists:

!----------------------------------------------------
! HISTORY (aka netCDF diagnostics)
!
! Loss of CO by OH for "tagged" species
!-----------------------------------------------------

! Units: [kg/s]
IF ( State_Diag%Archive_Loss ) THEN
   State_Diag%Loss(I,J,L,N) = Spc(N)%Conc(I,J,L) * k_Trop(2)     &
        * C(FixedOH) * DTCHEM
   !   C(ind_CO2_OH) / DTCHEM  &
   !   * State_Met%AIRVOL(I,J,L) * 1e+6_fp / XNUMOL_CO
ENDIF

How is k_Trop(2) derived and calculated? If I want to fix k_Trop(2), how should I modify it? My goal is to stabilize the GEOS-Chem CH4 lifetime. Should I modify a parameter in carbon_gases_mod.F90 or global_ch4_mod.F90? Could you give me some advice? Thanks.

yantosca commented 1 month ago

Thanks for writing @XunZhangNJU. The k_Trop(2) is the reaction rate in the KPP/carbon/carbon.eqn file as listed in this code snippet:

https://github.com/geoschem/geos-chem/blob/bef56c605e018eecbd91646a51ce82c7cd77f56a/KPP/carbon/carbon.eqn#L36-L59

k_Trop(2) is set in `KPP/carbon/carbon_Funcs.F90:

https://github.com/geoschem/geos-chem/blob/bef56c605e018eecbd91646a51ce82c7cd77f56a/KPP/carbon/carbon_Funcs.F90#L204-L207

from the function GC_OHCO which is located further down in KPP/carbon/carbon_Funcs.F90:

https://github.com/geoschem/geos-chem/blob/bef56c605e018eecbd91646a51ce82c7cd77f56a/KPP/carbon/carbon_Funcs.F90#L454-L482

Probably easiest way to set k_Trop(2) to a constant is to edit the GC_OHCO routine to be something like this:

FUNCTION GC_OHCO() RESULT( k )
    !
    ! Reaction rate for:
    !    OH + CO = HO2 + CO2 (cf. JPL 15-10)
    !
    ! For this reaction, these Arrhenius law terms evaluate to 1:
    !    (300/T)**b0 * EXP(c0/T)
    ! because b0 = c0 = 0.  Therefore we can skip computing these
    ! terms.  This avoids excess CPU cycles. (bmy, 12/18/20)
    !
    REAL(dp) :: klo1,   klo2,   khi1,  khi2
    REAL(dp) :: xyrat1, xyrat2, blog1, blog2,   fexp1
    REAL(dp) :: fexp2,  kco1,   kco2,  TEMP300, k
    !
    !---------------------------------------------------------------
    ! Hack to set reaction rate to a constant
    k = 1.000e-12_dp  ! NOTE: Replace this with the real reaction rate
    RETURN
    !---------------------------------------------------------------
    klo1   = 5.9E-33_dp * K300_OVER_TEMP
    khi1   = 1.1E-12_dp * K300_OVER_TEMP**(-1.3_dp)
    xyrat1 = klo1 * NUMDEN / khi1
    blog1  = LOG10( xyrat1 )
    fexp1  = 1.0_dp / ( 1.0_dp + blog1*blog1 )
    kco1   = klo1 * NUMDEN * 0.6_dp**fexp1 / ( 1.0_dp + xyrat1 )
    klo2   = 1.5E-13_dp
    khi2   = 2.1E+09_dp * K300_OVER_TEMP**(-6.1_dp)
    xyrat2 = klo2 * NUMDEN / khi2
    blog2  = LOG10( xyrat2 )
    fexp2  = 1.0_dp / ( 1.0_dp + blog2*blog2 )
    kco2   = klo2 * 0.6_dp**fexp2 / ( 1.0_dp + xyrat2 )
    k      = kco1 + kco2

  END FUNCTION GC_OHCO

XunZhangNJU commented 1 month ago

Thanks@yantosca. Actually, I want to set k_Trop(1) to a constant in the KPP/carbon/carbon.eqn file as listed in this code snippet:


       ! OH and Cl concentrations [molec/cm3]
       C(ind_FixedOH) = ConcOHmnd * OHdiurnalFac
       C(ind_FixedCl) = ConcClMnd

       ! CH4 + offline OH reaction rate [1/s]
       ! This is a pseudo-2nd order rate appropriate for CH4 + OH
       IF ( PCO_fr_CH4_use ) THEN
          k_Trop(1) = PCO_fr_CH4 * OHdiurnalFac
          k_Trop(1) = SafeDiv( k_Trop(1), C(ind_CH4)*C(ind_FixedOH), 0.0_dp )
       ELSE
          k_Trop(1) = 2.45E-12_dp * EXP( -1775.0E+0_dp /TEMP )  ! JPL 1997
       ENDIF
`
but I just couldn't get it right.

yantosca commented 1 month ago

Thanks for the feedback @XunZhangNJU. In that case you can do:

 // Rate k_Trop(1) [1/s] is set to HEMCO-supplied P(CO) from CH4 
 CH4 + FixedOH = CO    + COfromCH4   : 1.0e-12_dp  // replace with the actual rate law constant

Then rebuild the mechanism as described in

https://geos-chem.readthedocs.io/en/latest/geos-chem-shared-docs/supplemental-guides/using-kpp-with-gc.html

and then re-compile GEOS-Chem with

$ cd build
$ cmake -DMECH=carbon -DRUNDIR=..
$ make -j
$ make install

github-actions[bot] commented 3 days 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 issue from closing this issue.

geoschem / geos-chem