I need guidance about adding a species to the Carbon Simulation

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.

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.
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.
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?

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

Jared

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

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:

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].
If we aren't usingPCO_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

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.

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

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.

geoschem / geos-chem