I need guidance about adding a species to the Carbon Simulation

[JFBrewer]

Your name

Jared Brewer

Your affiliation

University of Minnesota

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

Hi Folks,

I'm trying to add ethane to the carbon simulation so that I can speed up my CHEEREIO workflows. As of now, I'm working in GC 14.4.2. To keep in line with what's currently in the Carbon sim, I'm only representing ethane loss to oxidation by OH and Cl (omitting Br and NO3, but these are considered minor loss processes). I'm having some problems getting this apparently simple simulation to work, however, and wanted to ask for some guidance from the experts.

WHAT I'VE DONE SO FAR:

• Added ethane emissions to HEMCO
• Confirmed that ethane emissions between the carbon sim and fullchem sim are similar
• Added the following two lines to carbon.eqn:
  C2H6 + FixedOH = C2H6P : 7.66d-12*EXP(-1020.0d0/TEMP);
C2H6 + FixedCl = C2H6P : 7.20d-11*EXP(-70.0d0/TEMP)*TROP;
• Added ethane-specific conversion factors to carbon_funcs.F90, specifically within carbon_ConvertKgToMolecCm3 and carbon_ConvertMolecCm3ToKg
• Added id_C2H6 and xnumol_C2H6 to the calls for the above functions in GeosCore/carbon_gases_mod.F90
• Added definition of xnumol_C2H6 and id_C2H6 and id_C2H6_adv to GeosCore/carbon_gases_mod.F90

Together, these fixes come close to getting the result I think I should see, but I'm still getting Ethane values that are much too high - 8.5x higher median ethane mixing ratios with the GCv5 ethane or 2.5x higher mixing ratios using OH from GCv14. The degree of overestimation is almost certainly too high to be simply the result of the missing oxidation terms from Bromine and NO3 - I think there still has to be a bug here.

I'm somewhat stumped for what other mistakes I might have made here, but I have a couple of ideas.

1. In the context of the carbon simulation, should I be thinking of ethane in terms of molecules ethane or molecules of carbon? I thought we'd moved away from the latter for the fullchem sim, but I realize I'm uncertain in the context of the carbon sim.
2. Do I need to make a change in how I'm formatting HEMCO lines when I add them into the Carbon sim? As of now, I'm just copying them over from fullchem, but it occurs to me that may not be correct.
3. Do I need to create another k_trop reaction for OH+C2H6, analagous to that for CH4+OH or CO+OH? If so, what should I use to represent the [OH]? I'm not 100% clear on the distinction between ind_FixedOH, ConcOHmnd, and OHdiurnalFac?
4. Is there any other set of custom code specific to the carbon sim that I've overlooked and need to add?

Thanks for your guidance - I'll keep plugging at this, too.

Jared

[yantosca]

Thanks for writing @JFBrewer.

In the context of the carbon simulation, should I be thinking of ethane in terms of molecules ethane or molecules of carbon? I thought we'd moved away from the latter for the fullchem sim, but I realize I'm uncertain in the context of the carbon sim.

The other species in the carbon simulation are in molec CH4, molec, CO, etc. You don't need to track them in molecules carbon.

Do I need to make a change in how I'm formatting HEMCO lines when I add them into the Carbon sim? As of now, I'm just copying them over from fullchem, but it occurs to me that may not be correct.

I think you can copy them over from fullchem. For example e.g. CEDS CO entries look like this:

https://github.com/geoschem/geos-chem/blob/bef56c605e018eecbd91646a51ce82c7cd77f56a/run/GCClassic/HEMCO_Config.rc.templates/HEMCO_Config.rc.carbon#L666C1-L681C1

which I think is more or less the same as in the fullchem.

Do I need to create another k_trop reaction for OH+C2H6, analagous to that for CH4+OH or CO+OH? If so, what should I use to represent the [OH]? I'm not 100% clear on the distinction between ind_FixedOH, ConcOHmnd, and OHdiurnalFac?
Is there any other set of custom code specific to the carbon sim that I've overlooked and need to add?

I think you can replicate what is done for CH4 + OH:

       !---------------------------------------------------------------------
       ! k_Trop(1): Rate [1/s] for rxn: CH4 + OH_E = CO + CO_CH4
       !---------------------------------------------------------------------

       ! 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

This is called from the routine Chem_Carbon_Gases:

       ! Compute the rate constants that will be used
       CALL carbon_ComputeRateConstants(                                     &
            I                = I,                                            &
            J                = J,                                            &
            L                = L,                                            &
            dtChem           = dtChem,                                       &
            ConcClMnd        = Global_Cl(I,J,L),                             &
            ConcOHmnd        = Global_OH(I,J,L),                             &
            LCH4_by_OH       = LCH4_by_OH(I,J,L),                            &
            LCO_in_Strat     = LCO_in_Strat(I,J,L),                          &
            OHdiurnalFac     = OHdiurnalFac(I,J),                            &
            PCO_in_Strat     = PCO_in_Strat(I,J,L),                          &
            PCO_fr_CH4_use   = Input_Opt%LPCO_CH4,                           &
            PCO_fr_CH4       = PCO_fr_CH4(I,J,L),                            &
            PCO_fr_NMVOC_use = Input_Opt%LPCO_NMVOC,                         &
            PCO_fr_NMVOC     = PCO_fr_NMVOC(I,J,L),                          &
            State_Met        = State_Met,                                    &
            State_Chm        = State_Chm                                    )

The ConcOHmnd and ConcClMnd are the monthly mean OH and Cl as read in from HEMCO. OH has a diurnal factor applied as well. If in your code the diurnal factor isn't being applied, that could potentially explain the overproduction.

Also tagging @msulprizio @lizziel

[JFBrewer]

Hi Bob - thanks for this response!

I remain just a bit confused about k_Trop, though. As I wrote in my first post, I currently set the C2H6 + OH rate in carbon.eqn, so it makes sense to me that I might be missing the OH diurnal scaling, which could explain the too-low C2H6 loss rate, in theory.

However, I'm confused about the idea of copying k_Trop for the CH4 + OH loss. As I read that code, it appears to me that there's two possible options for handling the CH4 loss:

1. If we've selected the PCO_fr_CH4_use option, then we read in some PCO files from the last benchmark via HEMCO and use that PCO rate x the OH diurnal factor divided by the product of [CH4] and [OH].
2. If we aren't using PCO_fr_CH4_use, then KCO is just an arrhenius expression, without any reference to the diurnal OH at all.

EDIT - sorry, I accidentally submitted too early. The fact that this interface and slack have opposite commands is a menace.

So how should I use this structure to add in the ethane reaction and take into account the diurnal OH term? Does it just look like case 2? Should there be some reference to the OH term in case 2/in my C2H6 reaction?

Thanks,
Jared

[JFBrewer]

To put a finer point on this, since I've just talked it through again:

I am interested in the possibility that I'm not applying the OH diurnal factor - that could make some sense as a source for my error - but where does the OH diurnal factor get applied? The only places I can find it used for Methane within carbon_Funcs.F90 is when the offline PCO_fr_CH4 option is used, and that is just to back out a K_CH4 term. If we don't use that, then we don't reference the OH Diurnal factor at all within carbon_ComputeRateConstants.

I'm trying to think about this in reference to my current strategy, which has been to use the Arrhenius expression within KPP, rather than K_trop. Does a call to K_Trop result in a different [OH] being used to calculate the effective-first-order rate coefficient than the use of an Arrhenius expression hard-coded into KPP? If not, should both my and the CH4 reaction rate when PCO_fr_CH4_use == FALSE have an OHdiurnalFac term appended to them, as in the case of k_Trop(3) (NMVOC + OH -> CO)?

Basically, I'm trying to figure out why my current strategy doesn't get the OH diurnal factor applied but a strategy analogous to the PCO_fr_CH4_use == FALSE case would.

[yantosca]

Hi @JFBrewer, thanks for the update.

but where does the OH diurnal factor get applied? The only places I can find it used for Methane within carbon_Funcs.F90 is when the offline PCO_fr_CH4 option is used, and that is just to back out a K_CH4 term. If we don't use that, then we don't reference the OH Diurnal factor at all within carbon_ComputeRateConstants.`.

The OH diurnal factor in the computing the rate stored CH4 + OH (k_Trop(1)), and FixedNMVOC + CO (k_Trop(3)).

Does a call to K_Trop result in a different [OH] being used to calculate the effective-first-order rate coefficient than the use of an Arrhenius expression hard-coded into KPP?

It shouldn't, at least that is my understanding. There is no OH used in the GCARR_* functions, these just take in the A, B, C constants for each reaction.

Tagging @msulprizio @msl3v

[msl3v]

Hi @JFBrewer
Could you try comparing the OH fields or even the rates (not rate constants) from the C2H6+OH reaction for a given gridbox between fullchem and your carbon sim version? My suspicion is that the OH is simply that much different between the offline and online sims.

[Liyl2]

Hi, @JFBrewer. I hope you’re doing well. I apologize for reaching out, but I have a similar need as you. I’m also looking to add some substances to the carbon mechanism. I’m currently running GC Classic 14.4.3, but I’m experiencing a hang during the carbon simulation (#2524 ). I successfully compiled the source code, but still can’t run the carbon mechanism properly.

Have you ever encountered this problem? I also tried GC Classic 14.4.2, while,but still same problem. I noticed that you’ve successfully run the carbon mechanism, and I was hoping you could take a look at my situation. I really appreciate your time! Here are the steps I followed:

[JFBrewer]

@Liyl2 It looks like there's some problem in your actual kpp. Did you successfully build the kpp mechanism with ./build_mechanism.sh after you made your changes to carbon.eqn and carbon_Funcs.F90? If you did, you could upload your carbon.eqn file here so I might have a better guess for you.

[Liyl2]

Thanks for your reply. I haven't made any changes to the carbon mechanism, and it won't work in the first step, just using the original source code. However, fullchem and CH4 simulation work well.

…

----------Reply to Message---------- On Thu, Oct 24, 2024 00:26 AM Jared ***@***.***> wrote: @Liyl2 It looks like there's some problem in your actual kpp. Did you successfully build the kpp mechanism with ./build_mechanism.sh after you made your changes to carbon.eqn and carbon_Funcs.F90? If you did, you could upload your carbon.eqn file here so I might have a better guess for you. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: ***@***.***>