Kinetic isotope effects of 12 CH 3 D + OH and 13 CH 3 D + OH from 278 to 313

Methane is the second most important long-lived greenhouse gas and plays a central role in the chemistry of the Earth’s atmosphere. Nonetheless there are significant uncertainties in its source budget. Analysis of the isotopic composition of atmospheric methane, including the doubly substituted species CH3D, offers new insight into the methane budget as the sources and sinks have distinct isotopic signatures. The most important sink of atmospheric methane is oxidation by OH in the troposphere, which accounts for around 84 % of all methane removal. Here we present experimentally derived methane+OH kinetic isotope effects and their temperature dependence over the range of 278 to 313 K for CH3D and CH3D; the latter is reported here for the first time. We find kCH4/kCH3D = 1.31±0.01 and kCH4/k13CH3D = 1.34±0.03 at room temperature, implying that the methane+OH kinetic isotope effect is multiplicative such that (kCH4/k13CH4)(kCH4/kCH3D)= kCH4/k13CH3D, within the experimental uncertainty, given the literature value of kCH4/k13CH4 = 1.0039± 0.0002. In addition, the kinetic isotope effects were characterized using transition state theory with tunneling corrections. Good agreement between the experimental, quantum chemical, and available literature values was obtained. Based on the results we conclude that the OH reaction (the main sink of methane) at steady state can produce an atmospheric clumped isotope signal (1(CH3D)= ln([CH4][ CH3D]/[ CH4][CH3D])) of 0.02± 0.02. This implies that the bulk tropospheric 1(CH3D) reflects the source signal with relatively small adjustment due to the sink signal (i.e., mainly OH oxidation).


Introduction
Atmospheric methane is the subject of increasing interest from both the climate research community and the public due its impacts on climate change, as reported by the IPCC (2013).The direct radiative forcing of methane is 0.64 W m −2 .Including feedback mechanisms and secondary effects (e.g., increased O 3 production, stratospheric water vapor, and production of CO 2 ), methane's radiative forcing becomes 0.97 W m −2 , two-thirds of the forcing by CO 2 over the same time period (IPCC, 2013, Fig. 8.15).
Atmospheric methane has both natural and anthropogenic sources and the two categories contribute about equally (Ciais et al., 2013, and references therein).Wetlands are the dominant natural source, and agriculture and waste are the largest anthropogenic sources.Fossil fuels make smaller contributions.The majority (84 %) of atmospheric methane is removed by oxidation by OH in the troposphere: The rest (4 %) is removed by soil (Kirschke et al., 2013).Carbon and hydrogen isotopic analyses are widely used to distinguish microbial and thermal sources of atmospheric methane (e.g., Lowe et al., 1997;Ferretti et al., 2005;Tyler et al., 2007;Lassey et al., 2007).However, Reactions (R1), (R2), and (R3) produce relatively large D/H isotope effects (Saueressig et al., 1995(Saueressig et al., , 1996(Saueressig et al., , 2001;;Crowley et al., 1999;Feilberg et al., 2005).Thus, the construction of an accurate top-down methane budget by isotopic analysis must take the isotopic signatures of both sources and sinks into account (Quay et al., 1999;Bergamaschi et al., 2000;Allan et al., 2001a, b).An isotope budget based on methane source (and sink) fractionations is an underdetermined systems (e.g., Pohlman et al., 2009).Recent advances in mass spectrometry (Eiler et al., 2013;Stolper et al., 2014) and laser infrared spectroscopy (Ono et al., 2014;Wang et al., 2015) facilitate measurement of rare doubly substituted isotopologues.The abundance of these "clumped" isotopologues (clumped refers to the rare isotopes being clumped together) generally follows a stochastic distribution (i.e., [ 12 CH 4 ][ 13 CH 3 D] = [ 13 CH 4 ][ 12 CH 3 D]).However, small deviations from stochastic distribution can be induced by thermodynamic (Ma et al., 2008;Stolper et al., 2014;Liu and Liu, 2016), kinetic (Joelsson et al., 2014;Wang et al., 2015), and photolytic processes (Schmidt et al., 2013;Schmidt and Johnson, 2015).Analysis of the clumped isotope anomaly in methane will yield unique constraints for the methane budget.Optical methods, as will be shown in this paper, provide high throughput and accuracy for overcoming the problems of analysis of clumped CH 4 .The difference and advantage of this approach is the additional information not available in single-isotope analysis, especially regarding the mechanism of formation, independent of the enrichment of D and 13 C in the starting material.
The kinetic isotope effect, j E α, is a characteristic property of each process: where i E is the most abundant (here, the lighter) isotopologue, j E the rare (heavy) isotopologue, and k(E + OH) the reaction rate coefficient for the reaction E + OH.As a measure of how much of a fractionation of 13 CH 3 D kinetic reactions produce, the apparent clumpiness, γ , is used.It is a measure of the effect of the clumped substitution on the reaction rate, as opposed to the combined effect of two single substitutions.It is defined as (Wang et al., 2015) γ A related measure is the ( 13 CH 3 D) value that quantifies the extent to which rare isotopes clump together to form a multiply substituted species, as opposed to a stochastic distribution (Ono et al., 2014): where and [ 13 CH 4 ] represent the concentrations of the different isotopologues.
The kinetic isotope effects for the singly substituted species CH 3 D and 13 CH 4 have been studied previously both experimentally and theoretically; see Tables 3 and 4, respectively.The kinetic isotope effect 13 C,D α for the reaction with OH is not described in the existing literature.The related kinetic isotope effect for the CH 4 + Cl reaction was measured at room temperature with the present setup by Joelsson et al. (2014) and found to be 1.60 ± 0.04.The kinetic isotope effects in the doubly substituted methane isotopologue CH 2 D 2 for the reaction CH 4 + Cl has been studied previously by Feilberg et al. (2005).
In the present study the kinetic isotope effects D α and 13 C,D α are determined using the relative rate method.Species concentrations in the reaction cell are determined using Fourier transform infrared (FTIR) spectroscopy.Further, D α, 13 C α, and 13 C,D α are calculated using quantum chemistry and transition state theory.

Experimental procedures
Sixteen experiments were conducted, numbered from 1 through 16 (see

Relative rate method
The experiments were carried out using the relative rate method on a semi-static gas mixture.The decaying concentrations of reactants were measured as a function of the extent of reaction.Considering two bimolecular reactions with second-order rate coefficients k A and k B , and assuming there were no other loss processes, the following relation holds: Here and [B] t represent the concentrations of compounds A and B at times 0 and t, respectively.The

Photoreactor
Experiments were carried out in the photochemical reactor at the University of Copenhagen, Department of Chemistry.The reactor has been described in detail elsewhere (Nilsson et al., 2009).It consists of a 100 L quartz cell with multipass optics surrounded by 16 UV-C fluorescent Hg lamps in a temperature-controlled housing.The cell is coupled to a Bruker IFS 66v/S FTIR spectrometer with either a mercury cadmium telluride (MCT) (Experiments 1-4) or an indium antimonide (InSb) detector (Experiments 5-16).The InSb detector has a better signal-to-noise ratio; the MCT detector is used in Experiments 1-4 for logistical reasons.Two thermocouple gauges are placed inside the temperaturecontrolled housing to monitor the temperature and a pressure gauge is connected to the cell to monitor pressure inside the cell.The temperatures and the pressures were logged every 0.5 s.

Laboratory procedure
Gas mixes were prepared by expanding H 2 O vapor (Milli-Q ultrapure water) into the chamber through a glass gas mani-fold.The two methane isotopologues CH 3 D (Experiments 1-8) (purity > 98 %, Cambridge Isotope Laboratories, Inc.) or 13 CH 3 D (Experiments 9-16), and CH 4 (purity > 99 %, Aldrich) and O 3 were flushed into the chamber with a N 2 buffer (purity 99.998 %, Air Liquide), all at the concentrations given in Table 1. 13 CH 3 D was synthesized using the Grignard reaction (see Joelsson et al., 2014).O 3 was generated from O 2 (purity 99.97 %, Air Liquide) using an ozone generator (model AC-20, O 3 Technology), preconcentrated before injection on silica gel cooled with ethanol and dry ice to −67 • C. The desired pressure in the cell (450 hPa) was obtained using N 2 as a bath gas.The starting pressure is chosen such that the pressure is high enough for the N 2 to quench O( 1 D) radicals, but low enough to keep the final pressure below atmospheric pressure.The gas mixture was left to rest for up to 1.5 h while several IR spectra were recorded to ensure that no instability or dark chemistry occurs in the gas mix.The UV-C lamps were lit for up to 5 min photolyzing at least 75 % of the O 3 according to O( 1 D) then subsequently reacts with H 2 O to yield OH: mainly due to O( 3 P) + O 3 .A pressure gradient was established and maintained throughout the filling process such that no gas leaked back from the chamber into the gas line.Spectra were recorded at each filling step.The procedure was repeated until the mix had a final pressure of 933 hPa.Two experiments were conducted at each of the temperatures 278, 288, 298, and 313 K for each of the two heavy methane isotopologues.Exact temperatures are listed in Table 1.After each experiment a dilution test was performed: 133 hPa was pumped out and the chamber is refilled with 133 hPa of N 2 .This was repeated 5-6 times.Ideally, concentration calculations from the spectral fits (data analysis described below) of the resulting spectra should give a linear fit with the slope of 1.The slope of these dilution tests is presented in Table 2.
In an extra experiment with 12 CH 3 D, N 2 O (Air Liquide, no purity information available) was added as an O( 1 D) tracer.
The results from this experiment are used as a benchmark to validate a model that was constructed to investigate the extent of O( 1 D) chemistry (see Sect. 2.5).An example of an experimental plot can be found in Fig. 4 and the full data set in Figs.S1-S8 in the Supplement.

Data analysis
The experimental IR spectra were analyzed using the program MALT, which simulates experimental FTIR spectra (Griffith, 1996) combined with nonlinear least-squares fitting to best-fit the calculated spectra to measured spectra (Griffith et al., 2012).The MALT program generates a simulated spectrum from initial estimates of the absorber concentrations and instrument parameters.The residual between the experimen- tal and simulated spectra is reduced through iteration.Simulated line shapes are generated using HITRAN absorption parameters (version 2008) (Rothman et al., 2009) convolved with an FTIR instrument function simulating the Bruker IFS 66v/S instrument.The InSb detector covers a spectral range from 1800 to 5000 cm −1 and the MCT detector covers a spectral range from 400 to 5000 cm −1 .The concentrations of 12 CH 3 D and 13 CH 3 D were calculated from spectral fits in the region 2140-2302 cm −1 (see Figs. 1 and 2).Interference from H 2 O, CO 2 , and CO was eliminated by including simulated spectra obtained from the HITRAN database in the fit.As there are no HITRAN data available for 13 CH 3 D in this region, the cross sections from 2000 to 2400 cm −1 for this isotopologue were estimated by shifting the spectrum of 12 CH 3 D (see Joelsson et al. (2014)).Concentrations of 12 CH 4 were calculated from spectral fits in the region 2838-2997 cm −1 .Interference from 13 CH 3 D was reduced by including temperature-adjusted reference spectra in the fit, and interference from 12 CH 3 D, H 2 O, and H 2 CO was by including simulated spectra obtained from the HITRAN database in the fit (see Fig. 3).The spectral windows were sometimes adjusted to exclude saturated lines.After [ x CH 3 D] (where x = 12 or x = 13) and [ 12 CH 4 ] were obtained from the spectral analysis, ln([ x CH 3 D] 0 /[ x CH 3 D] t ) was plotted against ln([ 12 CH 4 ] 0 /[ 12 CH 4 ] t ) as described in Sect.2.1.A straight line was fitted to these data points using a weighted total least-squares routine (York et al., 2004).The fitting procedure takes uncertainties in both dimensions into account.The uncertainty σ (ln([A] 0 /[A] t )) was calculated using standard error propagation: where σ ([A]) was obtained as output from MALT.The fitting procedure, performed using a MATLAB script (York et al., 2004), also yields an error estimation which is defined as the uncertainty of the kinetic isotope effect σ (α).An example of a straight line fit can be found in Fig. 4 and the full data set in Figs.S1-S8 in the Supplement.In the temperature dependence curve fitting procedure, the parameters A and B are from a linearized version of the Arrhenius equation: which is fitted to the experimental data.Also here, the method of York et al. (2004) was used.The temperature in the cell was taken as the spatial average of the measurements from two thermocouples inside the temperature housing.The experiment temperature was defined by the temporal mean of the spatially averaged temperature measurement series and the uncertainty of the experiment temperature was the standard deviation of the spatially averaged temperature measurement series.

Kinetic model
A kinetic model was used to determine the influence of O( 1 D), Reaction (R3), which rivals Reaction (R1).The model has been previously described and was used by Nils- The amount of CH 4 lost by Reaction (R3) can therefore be approximated by where k R8 = 1.27 × 10 −10 cm −3 s −1 and k R3 = 1.75 × 10 −10 cm −3 s −1 (Sander et al., 2011)

Theoretical procedure
Rate constants and kinetic isotope effects for CH 4 + OH were calculated using a procedure similar to that employed by Joelsson et al. (2014).

Computational chemistry calculations
The geometries of reactants, products, and transition states were determined using a geometry optimization procedure based on the unrestricted MP2 method (Møller and Plesset, 1934) and the aug-cc-pVQZ orbital basis set (Dunning Jr., 1989;Woon and Dunning Jr., 1993).Harmonic vibrational frequencies for all relevant isotopologues of reactants, products, and transition states were obtained at the same level of theory.The calculations were carried out using the Gaussian 09 program package (Frisch et al., 2009).
The electronic energy of the optimized structures were refined using the CCSD(T) method (Watts et al., 1993;Knowles et al., 1993Knowles et al., , 2000) ) with aug-cc-pVTZ, aug-cc-pVQZ, and aug-cc-pV5Z basis sets.The results from the different basis sets were used to extrapolate the electronic en-ergy at the complete basis set limit following the approach of Halkier et al. (1998) as described by Joelsson et al. (2014).

Rate constant calculations
The abstraction of a CH 4 hydrogen atom by OH can occur at four different sites.Depending on the methane isotopologue in question, these sites are either distinguishable or indistinguishable.Microscopic rate constants are calculated for hydrogen abstraction at each site using classical transition state theory with a tunneling correction factor, where E (i) is the reaction barrier height (including the zero-point vibrational energy) for the ith reaction path.η (i) tun is a tunneling correction factor obtained using the Wigner tunneling correction (Wigner, 1932).
TS and Q reac are the partition functions for the transition state and reacting pair, respectively.The total rate constant is obtained by summing over the microscopic rate constants.

Results and discussion
The results of the 16 individual experiments (eight for each isotopologue) are tabulated in Table 2.The resulting D α, 13 C α, and 13 C,D α values from the experimental and theoretical studies are tabulated in Tables 3, 4, and 5 along with previous experimental and theoretical results.The results are also shown in Figs. 5 and 6 for reactions of CH 3 D and 13 CH 3 D, respectively.The exponential curve fits yielded the parameters presented in Tables 3 and 5 giving D α Arr = 1.32 ± 0.13 and 13 C,D α Arr = 1.32±0.20 at T = 298 K.The mean of results of room temperature experiments Experiments 1 and 2 and Experiments 9 and 10 is D α exp = 1.31 ± 0.01 and 13 C,D α exp = 1.34 ± 0.03, respectively.It follows that γ exp = 1.02 ± 0.02 at T = 298 K, using 13 C α exp = 1.0039 ± 0.0002 (Saueressig et al., 2001), meaning that the clumped isotope might react slower relative to what would be predicted based on the kinetic isotope effects of CH 3 D and 13 CH 4 .However, if the Arrhenius parameters are used to calculate the kinetic isotope effects at T = 298 K, γ Arr = 1.00±0.18(i.e., the reaction has no clumping effect).All uncertainties are given as 1 standard deviation (σ ).The theoretical results give γ theory = 1.00 at T = 298 K.The present experimental room temperature results for D α agree, to within the error bars, with the previous experimental studies, with the exception of DeMore (1993).DeMore (1993) used FTIR spectroscopy with a slow flow setup where the two methane isotopologues were measured separately with a common reference compound.The low D α value in DeMore (1993) may be explained by interference from Table 3. Experimental and theoretical studies of D α.The temperature dependence studies are presented in the Arrhenius form D α(T ) = A(T /298 K) n exp(BT −1 ), where A, n, and B are tabulated; T is the given temperature range; and D α(T = 298 K) is the resulting kinetic isotope effect at T = 298 K.Where no B coefficient is presented, A can be taken as D α(T ) for the given temperature range.All uncertainties are given as 1 standard deviation (σ ).A relatively high rate of Reaction (R3) would reduce the final kinetic isotope effect, since the kinetic isotope effect for oxidation with O( 1 D) is smaller than the kinetic isotope effect for the oxidation of OH (Saueressig et al., 2001).The present experimental results of D α are also in good agreement with Masgrau et al. (2001) and Sellevåg et al. (2006).The theoretical calculations of D α give a slightly higher value than the experimental results, although they are in good agreement with the best Arrhenius curve fit at T = 298 K.The experimental room temperature result for 13 C,D α = 1.34 ± 0.03 agrees, to within the error bar, with both the theoretical value and the best estimate of the Arrhenius curve fit at T = 298 K.
The theoretical calculations show a very small temperature dependence; the variability in the experimental data is large compared to the value of the slope, making a quantification of the temperature dependence uncertain.Furthermore, the   theoretical analysis revealed that the primary cause for the kinetic isotope effect is the substantially reduced reactivity of the D atom, which, in turn, can be explained by a significant increase in reaction barrier due to changes in vibrational zero-point energy and to a lesser extent tunneling.

Atmospheric implications
At steady state, assuming no clumping in emissions, ( 13 CH 3 D) = ln(γ ).It follows that ( 13 CH 3 D) = 0.02 ± 0.02, implying that the clumped isotope effect of the OH reaction is small.Wang et al. (2015) studied a variety of natural samples and found ( 13 CH 3 D) values that vary from (−3.36 ± 1.42)‰ for samples of bubbling gases from The Cedars, Cazadero, California, USA, to (6.17 ± 0.34)‰ for samples from gas voids and hydrates in sediment cores from the northern Cascadia margin.Methanogenic bacteria in the laboratory ranged from −1.34 to 2.29 ‰.Natural gas was found to vary from 2.53 to 5.18 ‰.Given these natural variations in ( 13 CH 3 D) values of the sources, it is clear from isotopic mass balance that the fractionation from atmospheric oxidation by OH will likely have little effect on the observed ( 13 CH 3 D) value of tropospheric methane.Rather, this number will reflect the relative contributions of the sources.( 13 CH 3 D) would therefore be a more straightforward tracer for tracking methane sources than conventional isotopic analysis.Further field studies are clearly needed to characterize ( 13 CH 3 D) values of microbial and thermal methane.

Conclusions
We present experimentally derived CH 4 + OH kinetic isotope effects and their temperature dependence for CH 3 D and 13 CH 3 D; the latter is reported for the first time.We find D α = 1.31 ± 0.01 and 13 C,D α = 1.34 ± 0.03 at room temperature, implying that the kinetic isotope effect is multiplicative such that (k CH 4 /k13 CH 4 )(k CH 4 /k CH 3 D ) = k CH 4 /k13 CH 3 D to within the experimental uncertainty.We compare our experimental results to theoretical estimates derived using transition state theory with tunneling correction and kinetic isotope effects reported in the literature.We find good agreement between theoretical and literature values.Based on these experiments we find that the OH reaction (the main sink of methane) at steady state has a clumped ( 13 CH 3 D) = 0.02 ± 0.02.
The Supplement related to this article is available online at doi:10.5194/acp-16-4439-2016-supplement.

)Figure 1 .
Figure1.A typical spectral fit in the region where [ 13 CH 3 D] is obtained (Experiment 10).Experimental data are shown by the topmost line, followed by the fitted (synthetic) partial spectra of the most dominant absorbers ( 13 CH 3 D and H 2 O), the resulting fitted (synthetic) spectrum (also including CO 2 and CO), and the residual between the measured and fitted spectra.

Figure 2 .
Figure2.A typical spectral fit in the region where [ 12 CH 3 D] is obtained (Experiment 5).Experimental data are shown by the topmost line, followed by the fitted (synthetic) partial spectra of the most dominant absorbers (CH 3 D and H 2 O), the resulting fitted (synthetic) spectrum (also including CO 2 and CO), and the residual between the measured and fitted spectra.

Figure 3 .
Figure3.A typical spectral fit in the region where [ 12 CH 4 ] is obtained (Experiment 5).Experimental data are shown by the topmost line, followed by the fitted partial (synthetic) spectra of the most dominant absorbers (CH 4 and CH 3 D), the resulting fitted (synthetic) spectrum (also including H 2 O), and the residual between the measured and fitted spectra.

Figure 5 .
Figure 5.The individual measurements of D α are represented by black points, the accompanying individual error bars are represented by thin black solid lines, the experimental Arrhenius curve fit is represented by a black dashed line, the theoretical Arrhenius curve fit is represented by a blue dashed-dotted line, the mean of the room temperature measurements is represented by a red solid circle, the uncertainty of the room temperature measurements is represented by a thick red solid line, and the result from Gierczak et al. (1997) (which is included for comparison) is represented by a green dotted line.All uncertainties are given as 1 standard deviation (σ ).

Figure 6 .
Figure 6.The individual measurements of 13 C,D α are represented by black points, the accompanying individual error bars are represented by thin black solid lines, the experimental Arrhenius curve fit is represented by a black dashed line, the theoretical Arrhenius curve fit is represented by a blue dashed-dotted line, the mean of the room temperature measurements is represented by a red solid circle, and the uncertainty of the room temperature measurements is represented by a thick red solid line.All uncertainties are given as 1 standard deviation (σ ).

Table 1 .
Experimental setup.The experiment numbers are listed in column Exp., the detector is listed in column Detect., the heavy CH 4 isotopologue included in the experiments is listed in column [ x CH 3 D], the mean measured temperatures in the photoreactor are listed in column T , the H 2 O vapor concentrations at the start of the experiments (t = 0) as obtained from spectral fitting are listed in column [H 2 O] t=0 , the mean O 3 concentrations after refill (i.e., the "top" values) as obtained from spectral fitting are listed in column [O 3 ] top , the 12 CH 4 concentrations at the start of the experiments (t = 0) as obtained from spectral fitting are listed in column [ 12 CH 4 ] t=0 , and the heavy CH 4 concentrations at the start of the experiments (t = 0) as obtained from spectral fitting are listed in column [ x CH 3 D] t=0 .Note that, for the experiment including CH 3 D, the value of initial concentration only refers to [ 12 CH 3 D] t=0 .
t ) gives the relative reaction rate coefficient.In these experiments A is 12 CH 4 and B is 12 CH 3 D or 13 CH 3 D.

Table 2 .
Results.The experiment numbers are listed in column Exp., the heavy CH 4 isotopologue included in the experiments is listed in column x CH 3 D, the mean measured temperatures in the photoreactor are listed in column T , the kinetic isotope effect corresponding to the isotopologue is listed in column α, and the result of the dilution experiments is listed in column k dil .
*No dilution test performed.2.3 % [CH 4] lost by oxidation of O( 1 D), Reaction (R3).The kinetic model described above estimated that 4.7 % [CH 4 ] was lost by Reaction (R3) for this additional experiment.Both methods agree that O( 1 D) loss is a minor channel.No correction is applied, and the possible deviation is included in the estimated error.

Table 4 .
Experimental and theoretical studies of 13 C α.The temperature dependence studies are presented in the Arrhenius form13C α(T ) =A exp(BT −1 ), where A and B are tabulated, T is the given temperature range, and 13 C α(T = 298 K) is the resulting kinetic isotope effect atT = 298 K.Where no B coefficient is presented, A can be taken as 13 C α(T ) for the given temperature range.All uncertainties are given as 1 Tunable diode laser absorption spectroscopy-isotope ratio mass spectrometry.b Gas chromatography-mass spectrometry.c Isotope ratio mass spectrometry.d Transition state theory with Wigner tunneling correction.e Canonical variational transition state theory.f Conventional transition state theory.g Interpolated variational transition state theory with centrifugal-dominant, small-curvature tunneling coefficients.h Ab initio calculations.i No temperature information available. a

Table 5 .
Experimental and theoretical studies of 13 C,D α.The temperature dependencies are presented in the Arrhenius form 13 C,D α(T ) = A exp(BT −1 ), where A and B are tabulated, T is the given temperature range, and 13 C,D α(T = 298 K) is the resulting kinetic isotope effect at T = 298 K.All uncertainties are given as 1 standard deviation (σ ).Fourier transform infrared spectroscopy -relative rate.b Transition state theory with Wigner tunneling correction.c Average room temperature value, not obtained from curve fit. a