Supplement of Methane chemistry in a nutshell – the new submodels CH4 (v1.0) and TRSYNC (v1.0) in MESSy (v2.54.0)

Abstract. Climate projections including chemical feedbacks rely on state-of-the-art
chemistry–climate models (CCMs). Of particular importance is the role of
methane (CH4) for the budget of stratospheric water vapour (SWV),
which has an important climate impact. However, simulations with CCMs are, due
to the large number of involved chemical species, computationally demanding,
which limits the simulation of sensitivity studies. To allow for sensitivity studies and ensemble simulations with a reduced
demand for computational resources, we introduce a simplified approach to
simulate the core of methane chemistry in form of the new Modular Earth
Submodel System (MESSy) submodel CH4. It involves an atmospheric
chemistry mechanism reduced to the sink reactions of CH4 with predefined
fields of the hydroxyl radical (OH), excited oxygen (O(1D)), and
chlorine (Cl), as well as photolysis and the reaction products limited
to water vapour (H2O). This chemical production of H2O is
optionally fed back onto the specific humidity (q) of the connected
general circulation model (GCM), to account for the impact onto
SWV and its effect on radiation and stratospheric dynamics. The submodel CH4 is further capable of simulating the four most prevalent
CH4 isotopologues for carbon and hydrogen (CH4 and CH3D, as
well as 12CH4 and 13CH4). Furthermore, the
production of deuterated water vapour (HDO) is, similar to the
production of H2O in the CH4 oxidation, optionally passed back to
the isotopological hydrological cycle simulated by the submodel H2OISO,
using the newly developed auxiliary submodel TRSYNC. Moreover, the
simulation of a user-defined number of diagnostic CH4 age and emission
classes is possible, the output of which can be used for offline inverse optimization
techniques. The presented approach combines the most important chemical hydrological
feedback including the isotopic signatures with the advantages concerning the
computational simplicity of a GCM, in comparison to a full-featured
CCM.


In the simulation EMAC-apos-02, the CH4 submodel together with its isotopologue extension is applied. This includes isotopologues concerning both, carbon, and hydrogen isotopes. The submodel is set up with the KIEs as introduced in Table  1 (see main manuscript). The comprehensive interactive chemistry simulation EMAC-apos-03 is conducted with the kinetic chemistry tagging technique (MECCA_TAG) concerning hydrogen isotopologues, only. This configuration is chosen to investigate the pathways of deuterium from the source towards the end-product of deuterated methane (CH 3 D), i.e. deuterated water vapour (HDO). This requires to include KIEs for the intermediates, too, as well as to apply adequate branching ratios and isotope transfer probabilities. The inclusion of carbon isotopologues with MECCA_TAG is omitted due to the fact that MECCA_TAG introduces additionally nearly twice as many chemical reactions and species as included in the basic chemical mechanism. To maintain a computational efficient simulation, the CH4 submodel is in EMAC-apos-03 additionally applied to simulate the carbon related methane (CH 4 ) isotopologues. In this case, the CH 4 tracer of the simplified CH 4 chemistry (CH4) submodel (CH4_fx), acting as the master tracer for the CH 4 isotopologues in the CH4 submodel, is in each model time step reset to the CH 4 tracer in the Module Efficiently Calculating the Chemistry of the Atmosphere (MECCA) to ensure an identical overall CH 4 budget. The CH4 submodel also uses directly the on-line calculated the hydroxyl radical (OH), excited oxygen (O( 1 D)) and chlorine (Cl) distribution from MECCA.
The applied emission inventory in the presented simulations is an a posteriori inventory derived using an inverse optimization technique (Frank, 2018;Bruhwiler et al., 2005). The specific isotopic signatures of the emission sources used in the model are listed in Table S1.
The isotopic signatures are given in the δ-notation (McKinney et al., 1950). We use the standard isotopic signature of Vienna Standard Mean Ocean Water (VSMOW) for the signature of D in CH 4 (δD(CH 4 )) and Vienna PeeDee Belemnite (VPDB) for the signature of 13 C in CH 4 (δ 13 C(CH 4 )). Figure S1. Simulated multi-annual (2000Simulated multi-annual ( -2009 surface mixing ratio of CH4 in [nmoles of the chemical tracer per mole of air (mol mol −1 )] (upper), corresponding δ 13 C(CH4)VPDB in [‰] (middle), and δD(CH4)VSMOW in [‰] (lower). The left column shows results of EMAC-apos-02 and the right column those from EMAC-apos-03. The colored dots indicate the surface observations from NOAA/ESRL. The circles around the dots are the value of the simulation at the specific sampling height of the observation (in order to account for sub-grid orographic differences between simulation and observation).

S2.1 Surface sampling sites
To start with the evaluation of the simulation results, isotopic observations from NOAA/ESRL sampling sites (White et al., 2016(White et al., , 2017) are compared to the surface mixing ratios and δ-values of the simulations. For the comparison a climatological mean of 2000-2009 is used, since this time period is represented by most of the stations and the dynamic equilibrium of the simulated isotopic composition (as visible especially in EMAC-apos-03, Frank (2018)) has been reached.
EMAC-apos-03 agrees well with the stations regarding the CH 4 mixing ratio. Interesting is that the δ 13 C(CH 4 ) values are slightly better represented in EMAC-apos-02 compared to EMAC-apos-03, although the agreement is overall quite well in both simulations. This suggests that the emission signatures are a bit too low for methane containing 13 C ( 13 CH 4 ) in connection (a) (b) Figure S2. Taylor diagrams of the comparison between observations and the simulations EMAC-apos-02 (blue) and EMAC-apos-03 (purple) at various surface sampling sites. The Taylor diagram is shown for δ 13 C(CH4)VPDB (a) and for δD(CH4)VSMOW (b) with respect to the representation of the annual cycle during the considered time period [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009]. The size of the triangles indicates the bias in percent with upward oriented triangles indicating a positive and downward oriented triangles a negative bias, respectively. Circles indicate a bias of less than 0.1%. The symbols below the diagram are stations outside the displayed range of the Taylor diagram and are indicated by the colored number. The normalized standard deviation is displayed by the upper black number and the correlation coefficient by the lower black number on the right hand side of the symbol.
with the OH concentration in EMAC-apos-03. On the other hand, in case of δD(CH 4 ), EMAC-apos-03 agrees better, however, is still isotopically enriched compared to the station samples. This indicates that the chosen emission signatures for CH 3 D are too heavy.
In addition to that, the annual cycle of the observations is generally fairly well represented in both simulations (see Fig.  S2). However, the trend of the signatures at the stations over the years could not be captured yet. The reason for this is that the simulations fail to represent the general trend of the CH 4 mixing ratio and that the emission signatures of the individual sources are still uncertain.

S2.2 Airborne observations
During the CONTRAIL project, atmospheric air samples were taken with an Automatic air Sampling Equipment (ASE) mounted on a commercial aircraft (Umezawa et al., 2012). These air samples were later measured concerning the isotopic composition of CH 4 using a gas chromatography system and a flame ionization detector. The here presented sampling data comprise several flights between 2006 and 2010, with each flight providing up to 12 air samples.
The presented flights are seperated into two regions, as depicted in Fig. S3. The first region (green) indicates the flights on a north-south route, bound from Narita airport (Japan) to Sydney, Brisbane (Australia) or Guam, and the second region (red) represents those flights on an east-west route, bound from Narita to Honolulu (Hawaii).
Especially the first region provides the opportunity to investigate the representation of the meridional gradient and the northsouth imbalance in the δ-values in the model as it nicely spans over the tropics (40 • S−40 • N). Simulation results and the airborne observations in this region are depicted in Fig. S4, where green dots indicate the observations. The dark green line  Figure S3. Observations provided by the CONTRAIL project (Umezawa et al., 2012). The green shaded area indicates region 1, and the red shaded area indicates region 2.
(a) (b) (c) Figure S4. Comparison of airborne observations (green) in the meridionally aligned region 1 with simulation data from EMAC-apos-02 (blue) and EMAC-apos-03 (red). CH4 (a), δ 13 C(CH4)VPDB (b) and δD(CH4)VSMOW (c). The lighter red and blue colored markers indicate the de-biased simulation data for the direct comparison to the meridional gradient of the observations. The dark green line indicates the mean of the observations with the greenish shaded area being the corresponding single standard deviation.
indicates the mean of the observations and the shaded green area is the corresponding standard deviation. Simulated values are included as the red and blue dots respectively. It is apparent from the shown results that the meridional gradient in the simulations concerning CH 4 and both isotopic signatures are well represented, although the absolute values differ. This indicates that the implemented KIE in the model is reasonable and that adjustments to the signatures of the emission inventory are required.

S2.3 Balloon borne observations
The presented airborne observations are used to infer tropospheric chemical compositions. The high-altitude range of balloon borne observations enables to investigate the stratospheric isotopic signatures, as well.
The observational data are provided by Röckmann et al. (2011) and were obtained by altogether 13 balloon flights between 1987 and 2003 at four launch stations: Hyderabab in India (HYD), Aire sur l'Adour in France (ASA), Gap in France (GAP) and Kiruna in Sweden (KIR). The balloon-borne high-altitude air samples are obtained up to 10 hPa (35 km) and were later examined with respect to CH 4 mixing ratios as well as its isotopic composition concerning 13 CH 4 and CH 3 D using a high precision continuous flow isotope ratio mass spectrometer (Brass and Röckmann, 2010).
The observations shown in Fig. S5 indicate two features: -First, while CH 4 gets reduced towards higher altitudes, the isotopic content gets enriched (both, in δ 13 C(CH 4 ) and δD(CH 4 )). This occurs due to fractionation processes, which prefer lighter isotopologues in the sink reactions over heavier isotopologues.
-Secondly, again, a meridional gradient is visible. Polar regions tend to have less CH 4 than tropical regions, indicating to some extent the older age of the polar air masses. Consequently, polar regions are isotopically enriched compared to regions at mid and low latitudes.
The balloon-borne observations are compared to the simulations in Fig. S5 at pressure levels from 200 hPa to 10 hPa and separated into polar, mid-latitude and tropical regions. For the comparison, the monthly averaged data of the simulation is sampled at the specific year, month and location of the observation and interpolated from model levels to pressure levels. The plots in Fig. S5 further show the single standard deviation of the observations by the grey shaded areas and the standard deviation of all vertical profiles in the corresponding latitudinal region of the simulations as the shaded area in the color of the respective simulation.
The presented comparisons of observations to simulation results show that the global isotopic features of the meridional isotopic gradient and the isotopic gradient with altitude is captured well by both simulations (EMAC-apos-02, only with CH4 submodel, and EMAC-apos-03, with MECCA and the CH4 submodel). This indicates that the implementation of the simulation of CH 4 isotopologues is sufficiently realized and also confirms the suitability of the chosen KIE values. Absolute values and the inter-annual trend of the observations, however, are not captured well, which is mainly caused by uncertainties in the CH 4 emission fluxes and the applied source signatures. S3 Documentation of the CH4 submodel

S3.1 Introduction
The CH4 submodel represents a simplified CH 4 chemistry. It defines the tracer CH4_fx, which gets reduced via the four CH 4 sink reactions. The tracer is initialized from external data via the submodel TRACER (Jöckel et al., 2008) and modified by either emissions, which need to be introduced via the submodel OFFline EMISsions (OFFEMIS) (Kerkweg et al., 2006) or by Newtonian relaxation towards a lower boundary condition with the submodel TNUDGE (Kerkweg et al., 2006). Example namelist entries concerning the configuration of these submodels are found in Section S5. Additional to that, the CH4 submodel provides two further options. One is the simulation of the CH 4 isotopologues, and the second is the representation of age and emission classes of CH 4 , which, to some extent, are able to resolve an additional spatial and temporal information of the CH 4 emissions.
The option to simulate age and emission classes introduces additional tracers depending on the chosen number of age and emission classes. For every combination of age and emission class one tracer is defined, thus, if N is the number of age classes and M is the number of emission classes, in total N × M additional tracers are defined. The tracers are denoted by the names CH4_fx_e[mm]_a[nn], with [mm] being the identifying number of the emission class and [nn] the number of the age class.
The following section documents the subroutines, which are part of the CH4 submodel and in the section "User interface" the entries in the corresponding namelists are explained.

S3.2 MODULE messy_ch4_si: Subroutines in the submodel interface layer (SMIL)
These subroutines follow the general structure mandatory for Modular Earth Submodel System (MESSy) submodels. Note that _gp and _lg denote the Gaussian grid point and Lagrangian mode (see Brinkop and Jöckel (2019) for more information). In the presented examples solely the Gaussian grid point mode is used.
-SUBROUTINE ch4_initialize: Initializes the submodel, reads the control and coupling namelists and broadcasts the information to all parallel tasks.
-SUBROUTINE ch4_new_tracer: Defines the new tracers, which also includes the additional tracers regarding the submodel extensions (if applied).
-SUBROUTINE ch4_init_memory: Defines the channel objects and allocates memory.
-SUBROUTINE ch4_init_coupling: Sets pointers for coupling to the basemodel and other submodels.
-SUBROUTINE ch4_global_start: Sets values of internal variables with respect to the applied ageing method, if the option of age and emission classes is switched on.
S3.3 MODULE messy_ch4: Subroutines in the submodel core layer (SMCL) The following subroutines represent the core layer of the submodel.
proposed by Eichinger et al. (2015). The control (CTRL) namelist of the CH4 submodel includes the KIE values applied in the isotopologue extension of the submodel for all four sink reactions and both isotopologues. The KIE is represented in the form KIE = A·exp(B/T ), with A and B being the individual parameters and T the temperature in [K]. The namelist entries are given therefore as: KIE_CH4_XX_YY = A, B. XX and YY are set according to the specified reaction. XX denotes thereby the isotope in CH 4 and is 13C or D1. YY defines the reaction partner (either OH, O1D or CL) as well as the photolysis with jval. For those KIE, which are temperature independent, B is set to 0.0. The default values are A = 1.0 and B = 0.0, so that no KIE is applied.

S3.5.2 CH4 CPL namelist
The coupling (CPL) namelist of the CH4 submodel sets the parameters for the applied extensions and feedback on the specific humidity. It further determines the channel objects used as the reaction partners in the CH 4 oxidation. -L_GP and L_LG are both logical switches implying whether the Gaussian representation (GP) or Lagrangian representation (LG), or both are applied. The following namelist entries are shown for GP, however, there a identical entries for LG as well (indicated by gp and lg, respectively). (Default: L_GP = T, L_LG = F) -c_gp_OH, c_gp_O1D, c_gp_Cl and c_gp_jCH4 define the chosen channel objects for the reaction partners of CH 4 . They take two strings, the first indicates the channel, the second the object name.
-i_gp_nclass_emis_age denotes the number of emission-and age classes. It takes two integers, the first is the number of emission classes, the second is the number of age classes. (Default: i_gp_nclass_emis_age = 0, 0,) -r_gp_age_cll is an optional entry, which adjusts the time period (in days) of one age class. This entry is only valid for ageing option 1 and 2 (see main text section 3.1). (Default: 30.44 for each age class) -l_gp_adj_tend is a logical switch, which indicates whether the tendencies are adjusted so that the additional age and emission class tracers sum up to the master tracer CH4_fx. (Default: T) -i_gp_ageing is an integer switch indicating the ageing method, which means the advancing of CH 4 from one age class to the next older one. It can be chosen between: -0: monthly in one step

S4.1 Introduction
The submodel TRacer SYNChronization (TRSYNC) guarantees that the physical H 2 O tracers (incl. their isotopologues) receive also the correct tendencies of the corresponding chemical tracers.
The submodel CH4 defines the tracer HDO, the submodel H 2 O ISOtopologues (H2OISO) defines H2OISOHDOvap, and the MECCA_TAG in the MECCA defines I2H2O (or a different idiom, chosen by the user). The auxiliary submodel TRSYNC couples these tracers to combine the physical and chemical isotopic fractionation.
Without any isotopological extension solely the fifth-generation European Centre Hamburg general circulation model (ECHAM5) intrinsic tracer for specific humidity (q) is present. In this case, chemically produced H 2 O (either from CH4 or from MECCA) directly adds optionally to q. However, in case of an isotopological extension using H2OISO, CH4 and/or MECCA_TAG the following additional tracers are defined: -H2OISOHHOvap and H2OISOHDOvap (defined by H2OISO): The former is the total water tracer and the latter is the tracer of the rare isotopologue. Note that in H2OISO the two tracers do not add up to a master tracer, actually, H2OISOHHOvap represents and is identical to the master tracer (i.e. q).
-I1H2O and I2H2O, representing H 2 O and HDO, respectively (defined by MECCA_TAG): Both sum up to the chemical master tracer H2O.
-H2O (defined by MECCA): This tracer is originally not defined in MECCA, but is necessary in combination with MECCA_TAG for the internal scaling of I1H2O and I2H2O. Figure S6 depicts the schematics of the coupling. At the beginning of every time step, H2OISOHHOvap is set to the current value of q, correcting any numerical deviations of H2OISOHHOvap from q caused in the previous time step. Next, basically all tracers are modified by the same physical processes: advection, vertical diffusion and convection. However, for the submodels E5VDIFF, CONVECT and CLOUD the hydrological processes are doubled in H2OISO to allow for isotope effects. The submodel Multi-phase Stratospheric Box Model (MSBM) calculates a tendency for q, which is added to H2OISOHHOvap as well. An equivalent tendency is added to H2OISOHDOvap, which is derived such that no additional fractionation by the multi-phase stratospheric chemistry is implied.
After all physical processes are complete, the submodel TRSYNC is called. It takes care that all tendencies of the previous (physical) processes of HDO and I2H2O are deleted and overwritten by the corresponding tendencies of the H2OISO equivalent H2OISOHDOvap. I1H2O is exceptional, as it must be set to the difference of the total tracer H2OISOHHOvap and the rare isotopologue H2OISOHDOvap. Note that for technical reasons the tracer H2OISOHDOvap is defined as one half of the corresponding chemical isotopological tracers HDO and I2H2O.
Next CH4 computes the CH 4 oxidation and derives the feedback onto q and HDO. At the very beginning of MECCA, the intrinsic H2O tracer is synchronized with q. Before and after the calls of the kinetic solver, I1H2O and I2H2O are scaled appropriately to add up to H2O. After this, the feedback onto H 2 O is passed to q. To be precise, the sketch in Fig. S6 suggests that CH4 and MECCA are executed in the same simulation. This is indeed possible, but not necessary and it is important to note that only one of the two can provide the chemical feedback onto q, which can be arranged by corresponding switches in the namelists.
After the chemical processes, TRSYNC synchronizes the tracers HDO or I2H2O backward onto H2OISOHDOvap, and H2OISO also adds the chemical tendency of q to H2OISOHHOvap. As a last step H2OISO adjusts the tendency of H2OISOHHOvap so that it is conform to the tendency of q. The following section documents the subroutines, which are part of the TRSYNC submodel and in the section "User interface" the entries of the corresponding namelist are explained.

S4.2 MODULE messy_trsync_si: Subroutines in SMIL
These subroutines follow the general structure mandatory for MESSy submodels.
-SUBROUTINE trsync_initialize: Initializes the submodel, reads the coupling namelist and broadcasts necessary information to all parallel tasks.
-SUBROUTINE trsync_init_memory: Registers the tracers for the TENDENCY submodel, if the latter is applied.
-SUBROUTINE trsync_init_coupling: Sets pointers to the used tracers and checks whether the synchronized tracers are identical in terms of their molar mass.
-SUBROUTINE trsync_init_tracer: Initializes the tracers, hence checks whether the tracers are already initialized and accounts for a synchronized initial state.
-SUBROUTINE trsync_physc: This subroutine is called two times. The first time before the kinetic integrations of CH4 and MECCA and the second time after. It provides the necessary unit conversion and numerical adjustment to synchronize the chosen tracers.

S4.3 MODULE messy_trsync: Subroutines in SMCL
The following subroutines represent the core layer of the submodel. This subroutine converts the tendency tr_b_te from kg kg −1 moist air s −1 to mol mol −1 dry air s −1 .

S4.5.1 TRSYNC CPL namelist
The coupling (CPL) namelist of the TRSYNC submodel lists the tracers to be synchronized. TRSYNC takes two strings and one integer switch. The first string indicates the chemical tracer in mol mol −1 dry air . The second string indicates the physical tracer in kg kg −1 moist air . The integer string denotes, whether the synchronization is done in both ways (0), the chemical tracer is synchronized by the physical tracer before chemistry only (1), or the physical tracer is synchronized by the chemical tracer after chemistry (2).