of 29
Journal Pre-proofs
Theoretical and experimental constraints on hydrogen isotope equilibrium in
C1-C5 alkanes
Hao Xie, Michael J. Formolo, Alex L. Sessions, John M. Eiler
PII:
S0016-7037(24)00429-0
DOI:
https://doi.org/10.1016/j.gca.2024.08.023
Reference:
GCA 13537
To appear in:
Geochimica et Cosmochimica Acta
Received Date:
15 May 2024
Accepted Date:
29 August 2024
Please cite this article as: Xie, H., Formolo, M.J., Sessions, A.L., Eiler, J.M., Theoretical and experimental
constraints on hydrogen isotope equilibrium in C1-C5 alkanes,
Geochimica et Cosmochimica Acta
(2024), doi:
https://doi.org/10.1016/j.gca.2024.08.023
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover
page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version
will undergo additional copyediting, typesetting and review before it is published in its final form, but we are
providing this version to give early visibility of the article. Please note that, during the production process, errors
may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2024 Published by Elsevier Ltd.
1
Theoretical and Experimental Constraints on Hydrogen Isotope Equilibrium
in C1-C5 Alkanes
Hao Xie
a, b *
, Michael J. Formolo
c
, Alex L. Sessions
b
, John M. Eiler
b
a.
International Center for Isotope Effects Research, School of Earth Sciences and
Engineering, Nanjing University
b. Division of Geological and Planetary Sciences, California Institute of Technology
c. ExxonMobil Engineering and Technology Company
*Corresponding author. Email address:
hxie@nju.edu.cn
Abstract
Stable isotope ratios of C
1
–C
5
alkanes, the major constituents of subsurface gaseous
hydrocarbons, can provide valuable insights on their origins, transport, and fates. Equilibrium
isotope effects are fundamental to interpreting stable isotope signatures, as recognition of them in
natural materials indicates reversible processes and constrains the temperatures of equilibrated
systems. Hydrogen isotope equilibrium of C
1
–C
5
alkanes is of particular interest because
evidence shows that alkyl H can undergo isotopic exchange with coexisting compounds under
subsurface conditions. We present the results of a combined experimental and theoretical effort
to determine equilibrium hydrogen isotope distributions in mixtures of these hydrocarbon
compounds. We created two mixtures: one with C
1
, C
2
and C
3
(where C1 indicates methane, C
2
ethane, etc., hereinafter called ‘C
1
–C
3
mixture’) and another one with C
2
, C
3
,
i
C
4
,
n
C
4
,
i
C
5
and
n
C
5
(hereinafter called ‘C
2
–C
5
mixture’); in both cases, the mixtures were created to start out of
hydrogen isotope equilibrium prior to any subsequent experimental treatments. We tested the
performance of several metal catalysts as aids to H isotopic exchange by exposing these mixtures
to different metal catalysts at 100 or 200
̊C
and analyzing the compound-specific hydrogen
isotope ratios of the product gases. The C
1
–C
3
mixture exchanged hydrogen isotopes among the
starting gas components rapidly in the presence of Ru/Al
2
O
3
at 200°C. The isotope ratios reach a
steady-state (presumably equilibrium) after 72 hours of heating (up to 120 hours). The hydrogen
isotopic ratios of alkanes in the C
2
–C
5
mixture shifted significantly in the presence of Rh/Al
2
O
3
at 100°C and C
3
-C
5
compounds approached steady state after 70 hours, whereas C
2
had not yet
2
reached a steady state D/H ratio after 216 h of heating. We evaluated the reaction progress of
isotope exchange for each compound as a function of time with a reaction-network model
informed by constraints on chemical kinetics. The model indicates that C
3
,
i
C
4
,
n
C
4
,
i
C
5
and
n
C
5
in the C
2
-C
5
mixture reached internal isotope equilibrium after 70 h, hence we used their values
to calculate experimental equilibrium results. We also calculated equilibrium isotope
fractionations with the Bigeleisen-Mayer theorem using vibrational frequencies computed from
the density functional theory (B3LYP/aug-cc-pVTZ). We included torsional conformers and
explicit H positions in the calculations. We found that the experimental results from both the 200
̊C
C
1
–C
3
experiment and 100
̊C
C
2
–C
5
experiment are consistent with theoretical equilibrium
results within experimental uncertainty. Additionally, our reaction-network model for catalyzed
hydrogen isotope exchange between alkanes succeeded in fitting our experimental results. This
modeling framework can be adjusted to simulate exchange in natural settings for constraining
temperature histories and sources of natural gases.
Keywords: Alkane; Hydrogen Isotope; Natural gas; Thermodynamic Equilibrium
3
Introduction
C
1
–C
5
alkanes (methane, ethane, propane, butanes, and pentanes) in the natural
environment are important energy and chemical-engineering resources, potent pollutants, and
microbial metabolites. Stable isotope ratios of these alkanes have been useful in tracing their
sources, probing their formation pathways and environments, and discerning non-chemical
processes such as mixing and transport. In subsurface petroleum reservoirs, stable isotope ratios
of these compounds can indicate thermal maturity of gas generation (Berner and Faber, 1996),
type of source organic matter (Galimov, 2006), and microbial degradation or methanogenesis
(Whiticar, 1999; Boreham et al., 2008).
Most prior isotopic studies of the C
1
-C
5
alkanes have focused on carbon isotopes, with
relatively little attention focused on their hydrogen isotope ratios, especially for C
2
-C
5
species.
Far fewer samples have been measured for compound specific
2
H/
1
H than for
13
C/
12
C (Sherwood
et al., 2017; Milkov and Etiope, 2018) and less is known about the information that can be
deduced from the hydrogen isotope data that do exist. However, we value hydrogen isotope data
because it can track thermal evolution of alkanes by hydrogen exchange. Prior studies have
provided experimental and observational evidence that alkanes can be susceptible to catalyzed
H-isotope exchange inorganically at temperatures only slightly elevated compared to earth-
surface conditions (Sessions et al., 2004; Schimmelmann et al., 2006; Xie et al., 2020; Xie et al.,
2021) or by microbial metabolism (Turner et al., 2021). Hydrogen isotope exchange has the
potential to drive intermolecular (as well as intramolecular) isotope fractionations towards
equilibrium isotope distributions, which are controlled by thermodynamic equilibrium and
therefore suitable for quantitative geo-thermometry. Moreover, it has been shown that even non-
equilibrated isotopic distributions between and within hydrocarbons can usefully constrain
geological histories, as the progress of isotope exchange (or departure from equilibrium isotope
distributions) provides constraints on thermal maturity of natural gas (Xie et al., 2021).
Gauging the extent of isotope exchange reactions and attainment of equilibrium
distributions, and use of isotopic data from equilibrated samples for geothermometry all require
accurate constraints on equilibrium isotope effects (EIE; i.e., inter and intramolecular isotopic
distributions that are controlled by changes in free energies of molecules due to isotopic
substitution). More generally, equilibrium isotope effects are fundamental constants of broad
geochemical importance because they can be used to distinguish reversible processes from
4
irreversible processes. Therefore, knowledge of accurate equilibrium hydrogen isotope effects in
alkanes could serve as a cornerstone for interpreting hydrogen isotope data in studies not limited
to the topic of thermogenic natural gas, such as abiotic formation of hydrocarbons in
hydrothermal fluids and microbial hydrocarbon recycling in the surface aqueous environments
(Reeves et al., 2012).
EIE for hydrogen isotope substitutions in alkanes can be determined via two general
approaches: isotope analysis of materials that have been equilibrated under known, controlled
conditions, such as by hydrogen exchange reactions; and theoretical calculations using quantum
chemistry and statistical mechanics. We employed both methods in this study. The alkanes
examined in this study include the isomers of butane and pentane (excepting neopentane, which
is usually only a trace component of natural hydrocarbons and rarely analyzed for isotope ratios).
We conducted hydrogen exchange experiments on mixtures of C
1
–C
5
alkanes heated in
the presence of metal catalysts. On the surface of these catalysts, C–H bonds of adsorbed alkanes
can dissociate and recombine, allowing hydrogen to exchange across the various molecular
positions and metal surface (equation 1) (Bond, 2006; Sattler, 2018).
+
2푀⇌푅
+
(
1
)
In this reaction, R represents an alkyl group, M represents metal. This metal-catalyzed
exchange method has been previously shown to be effective for equilibrating hydrogen isotope
distributions in gaseous alkanes (Stolper et al., 2014; Xie et al., 2018). We analyzed the
compound-specific hydrogen isotope ratio analysis of product gases. Additionally, we calculated
the EIE using the Urey-Bigeleisen-Mayer theorem. Frequencies of fundamental vibrational
modes of the isotopomers of these alkanes are calculated with density functional theory (DFT).
Finally, we compared the experimentally observed and theoretically computed fractionations to
assess the self-consistency of experiment and theory, and as a basis for evaluating the accuracy
of theory as a constraint on isotopic properties and temperature conditions that have not yet been
observed experimentally.
Methods
Nomenclature
We report hydrogen isotope ratios as
δ
2
H values on the VSMOW scale. We report the
isotope fractionation between two compounds using the
ε
value:
5
2
=
2
1
2
1
(
2
)
Calculation of equilibrium isotope effects
Equilibrium isotope effects between gaseous alkanes are calculated with the Urey-
Bigeleisen-Mayer method (U-B-M; Bigeleisen and Mayer, 1947). Reduced partition function
ratios
values) are calculated for all of the singly-substituted and non-substituted isotopologues
of C
1
-C
5
alkanes from methane to
n
-pentane (excluding neopentane). The equilibrium isotope
effect between two alkanes is given by the ratio of their molecule-average
β
values, which are
here calculated as the average of the
β
values of the singly substituted isotopomers of each
compound. Harmonic vibrational frequencies used in U-B-M theory are computed using DFT
with the B3LYP hybrid functional and the Dunning correlation-consistent triple-zeta basis set
with added diffuse functions: aug-cc-pVTZ. We perform geometrical optimization and
vibrational frequency calculation at the B3LYP/aug-cc-pVTZ level using the Entos Qcore
quantum chemistry package (Manby et al., 2019). We use tight convergence criteria for self-
consistent field calculations by turning on the ‘tighter scf’ option in vibrational frequency
computation.
One of our previous studies (Thiagarajan et al., 2020b) reported calculated EIE of these
molecules using the same level of theory. However, that study only considered the most stable
torsional conformations of each alkane molecule, while certain alkane molecules (
n
-butane,
i
-
pentane, and
n
-pentane) possess multiple stable conformations. Moreover, that calculation used
one of the non-equivalent hydrogen positions on the same carbon atom as a representative value
for all hydrogens on that carbon. These two approximations might lead to errors; therefore, in
this study, we account for the temperature-dependent conformer distribution and non-equivalent
hydrogen positions. We created molecular geometry models for all stable conformers of
n
-
butane,
i
-pentane, and
n
-pentane (Figure S1). We calculate the equilibrium conformer
distributions and use them as weights for averaging the
β
values for each alkane molecule. The
fraction of a torsional conformers follows:
=
푅푇
푅푇
(
3
)
6
where
ΔG
i
stands for the Gibbs Free Energy difference between the i
th
conformer and a reference
conformer (the most stable conformer of each compound, specified in the Table S1), which is
approximated with the enthalpy and entropy values.
∆퐺
=
∆퐻
∆푆
(
4
)
ΔH
i
and
ΔS
i
, enthalpy and entropy differences, are also standardized to a reference
conformer. They are also calculated at the B3LYP/ aug-cc-pVTZ level. We report all of the
calculated
ΔH
i
,
ΔS
i
in Table S1. Vibrational frequencies for every isotopomer of interest are
available in the Supplementary Materials.
The U-B-M theorem used here involves several approximations that may influence our
results. The U-B-M theorem assumes harmonic vibration, rigid rotation, and the Born-
Oppenheimer approximation. Prior studies have shown that partition function ratios for hydrogen
isotopes can have nontrivial effects of vibrational anharmonicity and vibrational-rotational
coupling (Bigeleisen and Mayer, 1947; Richet et al., 1977). Corrections to the U-B-M theorem
and alternative
ab initio
methods have been developed to account for these effects (Liu et al.,
2010; Marklanda and Berneb, 2012; Webb and Miller, 2014; Webb et al., 2017). It has been
previously shown that the hydrogen isotope fractionation between methylene and methyl
moieties in propane, a system relevant to our research, suffers from a combined error ranging
from 1 to 3 ‰ at 0–200°C when calculated using the B-M equation (Webb and Miller, 2014; Liu
et al., 2021; Yin et al., 2024). Such errors are smaller than the long-term analytical uncertainty of
the method we use for compound specific
δ
2
H measurement (5‰), so we consider the U-B-M
theorem satisfactorily accurate for our study. Similarly, it was shown that harmonic theoretical
values are reasonably consistent (within 10‰ at 0–100°C) with experimental values for
hydrogen isotope equilibrium between water and ketone molecules (Wang et al., 2009a), in large
part because errors in
β
associated with the harmonic assumption tend to cancel out between the
two equilibrium partners.
Preparation of alkane mixtures in isotope disequilibrium
We created two alkane mixtures, one with methane, ethane and propane (hereinafter
called ‘C
1
–C
3
mixture’) and another with ethane, propane, iso-butane, n-butane, iso-pentane and
n-pentane (hereinafter called ‘C
2
–C
5
mixture’). The main reason for doing two sets of exchange
experiments separately is that hydrogen isotope exchange for ethane and higher alkanes have
7
different kinetics and interaction mechanisms with catalysts as compared to methane (Bond,
2006). On metal catalysts, methane usually exchanges more slowly than C
2+
alkanes (Sattler,
2018). Furthermore, our initial series of experiments for test catalysts and heating conditions had
revealed that C
4
and C
5
alkanes tend to be destroyed at 200°C (inferred from quickly diminishing
molecular abundances), which is likely a result of catalytic decomposition. As such, equilibration
of these compounds requires temperature to be lower than 200°C; however, under those
conditions exchange between methane and ethane in the presence of metal catalysts were too
sluggish. Consequently, we performed C
1
–C
3
experiments at 200°C and C
2
–C
5
experiments at
100°C.
Methane, ethane and propane were sourced from commercial pure gas tanks (99%) from
Air Liquide. Butanes and pentanes were from pure gas cylinders (99%) from Sigma Aldrich. All
components except methane were further cryogenically purified in a vacuum line at liquid
nitrogen temperature to remove air before mixing. We attempted to mix each component in an
amount that would contribute an equal number of hydrogen atoms to the overall mixture, as
opposed to mixing equal molar quantities (Table 1). The original methane had a
δ
2
H = -175‰,
which was too close to the expected equilibrium value in our C
1
-C
3
gas mixture. In order to
ensure that our experiments would provide unambiguous evidence of hydrogen isotope change in
methane, we added to this mixture a small amount of
12
C
2
H
1
H
3
(99%, from Sigma Aldrich)
spike. The final enriched mixture has
δ
2
H
C1
= -125‰. Both the C
1
–C
3
and C
2
–C
5
mixtures have
non-equilibrium
δ
2
H patterns exchange experiments (Figure 1).
Exchange experiments
Isotope exchange experiments were conducted by heating gas mixtures in the presence of
metal catalyst in Pyrex® tubes. Based on literature reviews of the activities of heterogeneous
catalysts (Bond, 2006; Sattler, 2018), we focused our study on five catalysts: Pd/C (10% loading,
powder), Pd/Al
2
O
3
(10% loading, powder), Rh/Al
2
O
3
(5% loading, powder), Rh/Al
2
O
3
(0.5%
loading, pellet), Ru/Al
2
O
3
(5% loading, powder). All catalysts are commercial products from
Sigma-Aldrich. Each tube (1–2cc) was loaded with 40-60mg of catalyst in an anaerobic
environment. Catalyst was then degassed by heating with a half-inch torch flame (500 – 600°C)
under vacuum for 3 minutes. One of the gas mixtures was then condensed into the tube at liquid
nitrogen temperature (-196
°C),
transferred through a glass vacuum line. We found that the metal
8
catalysts worked as an absorbent for trapping methane, which cannot be fully condensed by itself
at -196
°C.
We prepared 140 μmol of one or the other gas mixture in each sample tube. The tube
was then flame sealed and placed in an oven (Heratherm™ General Protocol Oven by Thermo
Scientific™) for a controlled time varying from several hours to weeks. The oven has automatic
temperature control that keeps the temperature within +/- 0.5 °C of the set point.
After heating, the Pyrex® tube was removed from the furnace and cooled to room
temperature in air. We then cracked the Pyrex® tube under vacuum using a flexible tube cracker
and transferred the evolved gas to a known volume in a glass vacuum line and recorded its total
pressure (constraining the amount of gas remaining after heating). The gas was then transferred
through the vacuum line to a second glass tube immersed in liquid nitrogen. For experiments
with C
1
–C
3
mixture, we put degassed 5A molecular sieve into the tube beforehand to trap
methane. This second glass tube was then flame-sealed and moved to the gas chromatograph -
mass spectrometer system (below) for isotopic analysis.
We conducted control experiments to rule out experimental artifacts, such as losses
and/or isotopic fractionations associated with gas handling or cracking or other side reactions. In
these experiments, mixed gases are loaded into the tube with catalysts but kept at room
temperature for 1 day.
Molecular and hydrogen isotope analysis
Molecular compositions of the mixtures before heating were analyzed at Caltech on a
quadruple gas chromatograph/mass spectrometer (GC/MS) system, the Thermo Fischer
Scientific ISQ. Chromatographic separation was accomplished using an Agilent GS-GasPro
capillary column (30m x 0.25mm ID x 0.25μm film thickness). Concentration of the C
1
-C
3
components are standardized using a reference gas mixture made of 80% methane, 10% ethane,
5% propane and 5% CO
2
. Concentrations of C
4
and C
5
alkanes are inferred from hydrogen
abundance based on m/z 2 peak area in the GC/pyrolysis/isotope ratio mass spectrometer
(GC/Py/IRMS) system.
Compound-specific hydrogen isotope analysis of starting and post-heating gas mixtures
was conducted at Caltech on a GC/Py/IRMS system, Thermo-Finnigan Trace GC Ultra and
Delta+XP IRMS with Conflo-3 interface. For this method, sample glass tube is cracked and gas
expanded into a small volume (~2ml) that has a rubber septum on one end. For samples of the
9
C
2
–C
5
mixture, this volume was warmed to 80 °C to avoid partial condensation of C
4
and C
5
alkanes, which might lead to isotopic fractionation of remaining vapor. For samples of the C
1
–C
3
mixture, the sample tube is heated at 150°C for two hours to desorb the alkanes from molecular
sieve (Stolper et al., 2014). We use a Hamilton gas-tight syringe to sample the gas from the
septum. Around 10–100μL of gas is used per injection. For measurements of the C
1
–C
3
mixture,
a GS-GasPro column is used. For measurements of the C
2
–C
5
mixture, either a GS-GasPro or a
Zebron ZB-5ms column is used. The pyrolysis tube used to convert alkanes to H
2
for isotopic
analysis is heated to 1375 °C. In-house H
2
gas or CH
4
gas injections bracket sample injections
and serve as reference standards. 1–3 replicate measurements are made for each sample. Long
term external precision for
δ
2
H is 5 ‰
(1σ
per injection).
Results
Computational results
The thermodynamic properties of the stable conformers of
n
-butane,
i
-pentane, and
n
-
pentane are listed in Table S1. We found that torsional conformers can exhibit small to moderate
EIE differences, up to 6 ‰ at room temperature (Figure 2). The weighted average of each alkane
is generally very close (less than 1 ‰) to the most stable conformer (t for
n
-butane, tg for
i
-
pentane and tt for
n
-pentane), which also has the highest abundance. Approximating the EIE of
full conformational configuration of an alkane with that of the most stable conformer is thus
accurate enough for C
1
–C
5
alkanes.
The EIE for non-equivalent H positions on the same carbon atom can differ more
substantially. In the g+g- conformer of
n
-pentane, the maximum EIE between two H positions on
the same carbon at 25 °C is 36 ‰ (Figure S2). In the tg
conformer of
i
-pentane, the maximum
equilibrium isotope fractionation between two H positions on the same carbon at 25 °C is 35 ‰.
Therefore, significant error could be introduced if it were assumed that all H positions on the
same carbon have the same EIE.
The molecular-average calculated equilibrium H isotope fractionations (standardized to
methane) is presented in Figure 3. At 300 K, the H isotope fractionation order from higher to
lower is
n
C
5
>
n
C
4
>
i
C
5
> C
3
>
i
C
4
> C
2
; at 500 K, the order is
n
C
5
>
n
C
4
>
i
C
5
> C
3
> C
2
>
i
C
4
.
Experimental results
10
In control experiments (where gas mixtures are handled as in full experiments, but no
heating occurs) using either Rh/Al
2
O
3
or Pd/C catalysts, we did not observe change in hydrogen
isotope compositions from those of the original mixtures (Table 2). These results support the
interpretation that isotope shifts in the heated experiments reflect isotope exchange (pontentially
with fractionations associated with alkane destruction through side reactions).
In experiments aimed at equilibrating the C
1
–C
3
mixture, we found that Ru/Al
2
O
3
is the
most active catalyst. At 200 °C, hydrogen isotope ratios changed significantly in a few days
(Figure 4). At the end of the longest of these experiments,
δ
2
H
C1
decreased by 30‰;
δ
2
H
C2
decreased by 27‰;
δ
2
H
C3
increased by 46‰. The hydrogen isotope compositions are not
distinguishable between the 72h sample and the 120h sample, indicating that the mixture had
attained a steady-state, and so plausibly reached equilibrium.
In experiments aimed at equilibrating the hydrogen isotope fractionations in the C
2
–C
5
mixture, we found that Rh/Al
2
O
3
(5% loading, powder) is the most active catalyst. At 100 °C,
hydrogen isotope compositions of every component shifted significantly in one day (Figure 5).
At the end of these experiments,
δ
2
H
C2
decreased by 28‰;
δ
2
H
C3
increased by 50‰;
δ
2
H
iC4
increased by 18‰;
δ
2
H
nC4
increased by 18‰;
δ
2
H
iC5
decreased by 19‰; and
δ
2
H
nC5
decreased by
6‰ (a negligible change, based on analytical uncertainties).
Δ
2
H values of C
3
,
i
C
4
,
n
C
4
,
i
C
5
and
n
C
5
alkanes values stabilized after the 77h sample, after which the range of
δ
2
H values was
within external analytical uncertainty (5‰) for each one of the C
3
–C
5
n
-alkanes.
Δ
2
H value of C
2
was on a steady decreasing trend throughout the experiment series, extending to the final sample.
Thus, these data are consistent with C
3
-C
5
alkanes having reached and maintained mutual
exchange equilibrium, but ethane did not (possibly due to slower exchange kinetics).
We report results of these two series of experiments in Table 2. In Figure S4 we report
results from several sets of experiments that did not undergo substantial hydrogen isotope
exchange. These sets of experiments either used lower temperature (100°C for C
1
–C
3
equilibration) or less effective catalysts (Pd/C).
Discussion
Computational results
We compare our results to a few relevant past studies that calculated EIE for alkane
molecules. Wang et al. (2009b) presented equilibrium hydrogen isotope effects for primary,
11
secondary and tertiary carbon positions based on U-B-M theorem with the B3LYP/6-311G**
level of theory on C
7
molecules, but their values might not be applicable to C
1
–C
5
alkanes
because local chemical structures are different. It has been discovered that Wang et al.’s values
are not suitable for butane isomers (Liu et al., 2021). It has also been demonstrated that accurate
‘cutoff’ calculations demand identical local structures within three bonds of the target position
(He et al., 2020). Therefore, EIE of positions in C
1
–C
5
alkanes cannot be estimated by results
from higher alkanes. Piasecki et al. (2016) reported EIE for isotopologues of methane, ethane
and propane calculated from the U-B-M theorem with the B3LYP hybrid functional and a split-
valence triple-zeta basis set: 6-311G**. Their results are almost identical to our results calculated
with the aug-cc-pVTZ basis set — e.g.,
ε
2
H
C3-C1
and
ε
2
H
C2-C1
values are within 0.2‰ between the
two studies at equivalent temperature.
Our theoretical results are potentially affected by the limitations of the U-B-M theorem
such as the higher order effects of vibrational anharmonicity and vibrational-rotational coupling.
Liu et al. (2021) calculated equilibrium isotope effects for isotopologues of n-butane and i-
butane using the corrected U-B-M theorem with the coupled cluster method, CCSD(T), and the
6-311G** basis set. They applied a series of corrections to the U-B-M theorem, in order to
account for vibrational anharmonicity, vibrational-rotational coupling, quantum-mechanical
rotation, centrifugal distortion and hindered internal rotation and diagonal Born-Oppenheimer
correction. For
ε
2
H between n-butane and i-butane, our results are in good agreement (<0.4 ‰ at
100 °C; Table 3) with results from uncorrected and corrected CCSD(T)/6-311G** values.
Additionally, we can compare our results of hydrogen isotope fractionation between the
molecular positions of propane with prior studies that employed path-integral Monte Carlo
(PIMC) or the corrected U-B-M theorem (Yin et al., 2024). The central–terminal fractionation
factor of propane in our calculations differs by < 5 ‰ from both prior studies at 100 °C.
Therefore, the magnitude of higher order inaccuracies of the U-B-M theorem is smaller than
analytical error for the relevant molecules. However, we note for the purposes of possible future
studies that high-resolution mass spectrometry is capable of measuring
2
H of alkanes with
precisions as good as ~0.1 ‰, and so it is possible that one could experimentally test the relative
accuracies of these various models.
In an earlier work (Thiagarajan et al., 2020b), we calculated EIE of these molecules using
the same level of theory (B3LYP/aug-cc-pVTZ) as our study but only accounted for the most
12
stable torsional conformers and only one H atoms per C position. We found our updated results
can differ from that study (Figure S3). In
i
C
4
, the difference can be as high as 10‰ at room
temperature. As noted above, this change is caused by accounting for EIE variations in H
positions on the same carbon.
Experimental isotope fractionation values at steady state
In the C
1
–C
3
exchange experiments, the final sample at 120 h does not differ in isotope
composition from the 72 h sample, suggesting that steady-state, and therefore plausibly
equilibrium isotopic fractionations has been reached. As such, equilibrium isotope effects among
these three species can be inferred from the mean of the equilibrated samples of the time series.
The equilibrium isotope fractionation values from this experiment,
ε
2
H
C1-C3
and
ε
2
H
C2-C3
, are
consistent with the computational results (Figure 6). We standardize all
ε
values in the exchange
experiments to propane, because propane is present in both series of experiments.
The C
2
–C
5
experiment showed that most of the compounds reached steady isotope ratio
after 48 hours of exchange, except for ethane, which is still trending downwards in isotope ratios
of the last samples. We conclude that the subgroup, C
3
,
i
C
4
,
n
C
4
,
i
C
5
and
n
C
5
, reached mutual H
isotope steady state relatively quickly. This is perhaps due to faster exchange kinetics of these
higher alkanes than ethane rather than side reactions of ethane because the molecular
composition remained relatively stable during the course of this experiment (Table S2).
However, the coexistence of this subset with ethane that continues to evolve in hydrogen isotope
ratio implies that the system as a whole never reached full internal isotopic equilibrium. This
raises the question of whether the compounds that reached steady state had equilibrated with
each other, or instead have compositions that reflect a quasi-steady-state that differs from
equilibrium. The latter scenario may result from differential exchange of molecular positions,
where some positions equilibrate quickly while slower positions remain dis-equilibrated. We
proceed by estimating the time threshold for attainment of equilibrium among the more
exchangeable subset of compounds, by describing the progress of exchange reactions using a
chemical kinetic model that tracks every molecular position.
Modeling isotopic evolution during the hydrogen isotope exchange process
13
A common mathematical treatment of isotopic exchange is a pseudo-first order
approximation, which has isotopic ratios or difference between isotopic ratios approach the
equilibrium value via an exponential decay function (Roberts, 1939; Criss et al., 1987; Sessions
et al., 2004; Labidi et al., 2020). However, this simplistic form does not apply to a multi-
endmember exchange environment like our experiments. In our experiments,
δ
values and
ε
values can have more complicated paths to equilibrium because different compounds (and
moieties within each compound) undergo hydrogen isotope exchange at different rates. The
variability in exchange rates can be partially attributed to a decrease in bond dissociation
energies (BDE) with increasing carbon degrees (Sattler, 2018).
Therefore, we develop a more sophisticated chemical kinetic model to simulate the
evolution of hydrogen isotope compositions in our experiments. This model attempts to build the
full reaction network of catalytic exchange, such that we can evaluate the progress of hydrogen
isotope exchange and constrain the equilibrium values that would be attained at longer times.
This model considers the elementary reactions of adsorption (forward) and desorption
(backward) of alkanes on a metal surface (Equation 1). For example,
1
3
2
+
2푀⇌
1
3
+
2
(
5
)
Note that this formulation represents the uni-adsorption mechanism. Alternative
mechanisms such as multi-adsorption exists for C
2
and higher order alkanes, so this treatment is
an approximation. Multi-adsorption refers to the case where multiple carbon positions within a
single molecule are simultaneously bonded to metal atoms. This can occur in various
configurations, such as
α-β
multi-adsorption, where two adjacent carbon atoms are bonded;
α-γ
multi-adsorption, where the bonding involves two carbon atoms separated by one other carbon
atom (i.e., two bonds apart); and
α-β-γ
multi-adsorption, where three consecutive carbon atoms
in the molecule are bonded to metal atoms. For example,
α-β
adsorption is observed as a main
mechanism for ethane exchange on platinum surface (Loaiza et al., 1996). Unfortunately, there
has been no constraints of the relative strengths of these mechanisms for the experimental
conditions (catalysts, temperature) that we used, so we do not include them in our model.
The molecules our model considers include gas phase molecules with no isotopic
substitutions and single-deuterium substitutions, metal-bound molecules with no substitutions,
and metal bound
1
H and
2
H. The reactions include forward and backward versions of eqn 3,
which is automatically generated based on a list of isotopologues (included in the Supplementary
14
Materials). This model has a total of 44 unique isotopologues of the various molecular species
and 60 unique isotope exchange reactions. The kinetic treatments of reactions are based on the
observation that the activation energy of the desorption reaction is proportional to the BDE of
that C–H bond. We use methane adsorption/desorption as a reference reaction, to normalize
kinetic rate constants of reactions for other molecules. The rate constants of a non-methane
reaction could then be calculated from
0
=
exp
푠푓
훥퐵퐷퐸
푅푇
(
6
)
where
푠푓
denotes a scale factor, k
0
denotes rate constants for methane adsorption, R is ideal gas
constant and T is temperature in Kelvin. BDE is 438.9 kJ/mol for methane, 427.3 kJ/mol for
ethane and 420.9, 410.5, and 400.0 kJ/mol or primary, secondary and tertiary carbon positions of
C
3
to C
5
alkanes, respectively (Gribov et al., 2003; Luo, 2007).
For the
2
H-substituted version of reaction equation (1), we modify the reaction rate
constant with a kinetic isotope effect (KIE). We estimate the KIE =
exp(-ΔE
a
/RT), where the
activation energy difference
(ΔE
a
) is approximated as the zero-point energy (ZPE) differences
(ZPE
1H
-ZPE
2H
) in the reactants, as we are not aware of a previous quantitative simulation on the
isotope substitution in the transition states of adsorption. In practice, we set KIE (k
2
H/k
1
H) to
equal to
1/β,
where
β
is the reduced partition function ratio. This treatment also directly
parameterizes the equilibrium isotope effects into the kinetic model.
We initialize the model with all alkane molecules in the gas phase (adsorbed molecules
and hydrogen atoms are set to have zero concentrations). The starting isotopic compositions of
the gas molecules equals those of the initial gases in the prepared mixtures (Table 1). Initial
intramolecular (position-specific) isotope distributions are assumed to be homogenous except for
propane, which is assigned with
δ
2
H
Center
= -204 ‰ and
δ
2
H
Terminal
= -171 ‰ as they were
determined in a previous study (Xie et al., 2018). The reaction network is modelled via a set of
ordinary differential equations. We solve it numerically with a variable-step, variable-order
solver based on the numerical differentiation formulas of orders 1 to 5 (MATLAB ode15s).
We fit parameters of the reaction network model to experimental data by minimizing the
sum of squared error using the @fmincon function in MATLAB. We used the ‘interior-point’
algorithm for the C
1
–C
3
experiment and the ‘active-set’ algorithm for the C
2
–C
5
experiment. The
error is calculated as difference between the
ε
values from the model and those from the
15
experiments. We fit
2
H/
1
H
β
values of the first position (IUPAC naming rules) of each compound
in the exchange experiments. For the other positions within the same compound, we calculate
2
H/
1
H
β
values by combining
2
H/
1
H
β
values of the first positions with theoretical
β
ratios
between positions (provided the Supplementary Material). Additionally, we fit the scale factor
for converting BDE to activation energy (as in equation 5) as well as the rate constant of the
reference reaction
k
0
, which is the methane adsorption reaction. The code and data of this
exchange kinetic model are available at
https://github.com/1995123xh/alkexchange.
The exchange kinetic model succeeded in fitting our experimental data for both
experiments (Figure 7). In the C
1
–C
3
experiment, the reduced chi-squared is 32.8; in the C
2
–C
5
experiment, the reduced chi-squared is 30. The high goodness of fit between model trends and
experimental data on the non-equilibrated part of the time series supports our assumption that
exchange kinetics can be approximated by a function of BDE of C–H bonds. This simplification
assumes C–H bond dissociation is the rate limiting step of all exchange steps, which is consistent
with the most commonly accepted
Horiuti−Polanyi
mechanism of catalytic hydrogenation and
dehydrogenation (Sattler et al., 2014). We note that various other physicochemical phenomena
might also have influence on alkane exchange rates. Such phenomena might include electron
density distribution (Sattler et al., 2020), formation of
π-bonded
alkene intermediates (Natal-
Santiago et al., 1997), strength of physical adsorption and stereochemical hinderance (Clarke and
Rooney, 1976). The study of these higher order catalytic effects is beyond the scope of our work.
The exchange model showed that most of the alkane fractionation values follow a
monotonic trajectory (i.e. distance to equilibrium decreases with time). However, the
n
C
4
–C
3
fractionation displayed slight decrease at the beginning that is followed by an upward trend. This
might reflect the expression of different exchange kinetics of –CH
3
and –CH
2
– groups in
n
C
4
and
C
3
. The
ε
values each alkane pairs in the model (except for C
2
in the C
3
–C
5
experiment, which
continued to evolve) reached a stable steady state, which aligns with the equilibrium
fractionation factors embedded in the model. Therefore, the exchange model indicates that the
steady state observed in our experiments can be considered as isotopic equilibrium. In the 200
°C
C
1
–C
3
experiment series,
ε
values from the model approached equilibrium after 60 h, thus the
120 h sample and the 72 h sample should be indistinguishable from equilibrium. In the 100
°C
C
2
–C
5
experiment,
ε
values from the model approached equilibrium after 70 h, thus samples
between 77–216 h should be indistinguishable from equilibrium. Based on this result, the
16
experimental equilibrium fractionation factors can be calculated as the mean
ε
values for the
equilibrated samples in each of the 200
°C
and 100
°C
series respectively.
We compare our experimental results on equilibrium isotope effects with theoretical
predictions in Figure 6. They are overall in good agreement (smaller than or near
of
experimental results) with B3LYP results for both series of experiments (100°C and 200°C). For
the molecule pair of
n
C
4
i
C
4
, the theoretical value of
from Liu et al.’s calculation using
CCSD(T)/6-311G** with B-M corrections (Liu et al., 2021) is also consistent with uncertainty of
the experimental result. This agreement demonstrates that results from U-B-M theorem with
harmonic vibration and rigid rotation are adequate given the current precision of online
GC/Py/IRMS analysis. Alternative analytical techniques can reduce measurement uncertainty —
the offline reduction method (Schimmelmann, 1991) can achieve an uncertainty of 1–2‰ and the
direct molecular method (Stolper et al., 2014) can achieve an uncertainty of 0.1–0.2‰. When
these more precise techniques are applied, the nuanced differences in theory will be important.
Implications
Prior analysis of compound-specific hydrogen isotope compositions of C
1
-C
3
compounds
in natural gases have shown that they are generally not equilibrated at lower thermal maturity
and equilibrated at higher maturity (Xie et al., 2021). Therefore, the distance from intermolecular
hydrogen isotope equilibrium can potentially serve as a maturity indicator. The distance is
calculated as
|
푒푞푢푖푙
|
, where
denotes hydrogen isotope compositions of a gas sample (e.g.
2
H
C1
,
δ
2
H
C2
,
δ
2
H
C3
>, and
푒푞푢푖푙
represent that of this gas sample after full equilibration.
푒푞푢푖푙
value can be calculated using EIE obtained from this study. As thermal maturity increases, this
distance diminish and approach 0 at the dry gas stage.
This study demonstrates that theoretically calculated equilibrium isotope effects are
consistent with experimental ones, strengthening the foundation of natural sample interpretation.
When calculating hydrogen EIE using quantum chemical methods, our results show potentially
large difference between hydrogen atoms connected to the same carbon but have different C–H
orientations. Therefore, it is essential to take into account the positions of all hydrogen atoms
across the various torsional conformers. Otherwise, errors can emerge.
Additionally, the same argument now can be expanded to consider the evolution in
hydrogen isotope compositions of slightly larger molecules of
i
C
4
,
n
C
4
,
i
C
5
and
n
C
5
, which are
17
largely absent in dry gases but present in wet gas and condensate. Both theory and experiment
suggest that there are sizable equilibrium
2
H/
1
H fractionations between the
normal-
and
iso
-
isomers of the same carbon number (i.e.
n
C
4
/
i
C
4
or
n
C
5
/
i
C
5
) (Figure 6). The magnitude of
equilibrium isotope fractionation between isomers can be larger than kinetic fractionation values
generated by cracking reactions, because the magnitude of (molecular average) kinetic
fractionation associated with cracking is likely controlled by the number of carbons in the
product molecule (Chung et al., 1988). The discrepancy between equilibrium and kinetic isotope
fractionations in these isomer pairs (i.e.
n
C
4
/
i
C
4
or
n
C
5
/
i
C
5
) thus allows determination of
equilibrium and non-equilibrium controls, which can be used for tracking thermal maturation as
previously described.
The success of our isotope exchange kinetic model based on elementary reactions in
describing the trajectory of hydrogen isotope ratios in the exchange experiments supports the
validity of this model’s key assumption, that activation energy for hydrogen isotope exchange is
proportional to the BDE of C-H bonds. This model also shows that certain
δ
2
H and
ε
2
H values do
not approach the equilibrium value in a straightforward, monotonic way at all times. This
modeling framework can be applied to study exchange behavior in natural systems. We
exemplify the application of this exchange model by simulating the hydrogen isotope trajectory
starting from an immature natural gas sample (from the Julia gas field in the Walker Ridge Basin
(Thiagarajan et al., 2020a)) towards equilibrium (Figure 8; data and code available at:
https://github.com/1995123xh/alkexchange). Such analysis can be used to constrain the
temperature-time history of hydrocarbon fluids, or to determine the pristine cracked gas
characteristics before exchange takes place.
However, users should be aware of the limitations on the utility of this exchange kinetic
model in natural settings. First, hydrogen isotopes of gas compounds can reflect sustained
cracking or degradation rather than (or in addition to) hydrogen exchange reactions after
generation. This issue can be addressed by integrating kinetic isotope fractionation of
thermogenic gas generation and hydrogen exchange into the full model (Xia and Gao, 2019; Xie
et al., 2021). Additionally, we note that hydrogen isotope exchange can be catalyzed by
mechanisms other than metal-organic bonding, and the BDE based rate estimation used in our
model might not apply to other scenarios. For example, catalysis by metal oxides (Al
2
O
3
) is
suggested to take the carbanion mechanism that favors exchange on lower degree carbons (i.e.
18
acidity controls) (Sattler, 2018), which is opposite from BDE control. Catalysis by clay minerals
is suggested to take the carbocation mechanism, which favors exchange on the position adjacent
to carbocation-stable positions (Alexander et al., 1984). However, the H-abstraction reactions
have the same kinetic trend (tertiary>secondary>primary) as BDE control (Tsang, 1990). Such
differences in catalytic mechanisms will need to be considered when modeling natural systems.
Nevertheless, the modeling framework we present can be adjusted to reflect changes in kinetic
rules for these scenarios.
Data Availability
Data are available through Mendeley Data at: https://data.mendeley.com/datasets/bxp6k4cgf9/1
CRediT authorship contribution statement
Hao Xie:
Conceptualization, Methodology, Validation, Formal analysis, Investigation, Data
curation, Writing – original draft, Writing – review & editing, Visualization, Project
administration.
Michael J. Formolo:
Resources, Writing – review & editing, Funding
acquisition.
Alex L. Sessions:
Resources, Writing – review & editing, Supervision.
John M.
Eiler:
Conceptualization, Resources, Writing – review & editing, Supervision, Funding
acquisition.
Acknowledgements
This work was supported by Agouron Geobiology Postdoctoral Fellowship, Caltech,
ExxonMobil Upstream Integrated Solutions, the Caltech Joint Industry Partnership for Petroleum
Geochemistry, and the Department of Energy BES program. We thank Camilo Ponton, Fenfang
Wu, Nami Kitchen and Guannan Dong for advice on experiments and measurements. We thank
Brian Peterson for discussion of quantum chemical calculations.
Figure Captions
Figure 1. Initial hydrogen isotope compositions of our hydrocarbon mixtures. Left: C1–C3
mixture; right: C2–C5 mixture.
19
Figure 2. Equilibrium H isotope fractionation between the torsional conformers of
n
-butane,
i
-
pentane and
n
-pentane, normalized to the weighted average of each compound at conformational
equilibrium. g: gauche; t:trans; gmgm: gauche- gauche-; gmgp: gauche- gauche+.
Figure 3. Equilibrium H isotope fractionation between the
n
-alkane compounds and methane.
Figure 4. Compound-specific hydrogen isotope compositions in the C1–C3 experiment at
200
°C,
using Ru/Al
2
O
3
as a catalyst.
Figure 5. Compound-specific hydrogen isotope compositions in the C2–C5 experiment at
100
°C,
using Rh/Al
2
O
3
as a catalyst. Left: C3–
n
C5 alkanes; right: C2 ethane.
Figure 6. Comparing experimentally determined equilibrium hydrogen isotope fractionation for
alkane pairs with theoretical values. Each pair is normalized to propane (C3), except for the pair
of
n
C4–
i
C4, which is presented to compare with the anharmonic result (CCSD(T) corr) in Liu et
al. (2021).
Figure 7. Hydrogen isotope fractionation factors of alkanes in the exchange experiments,
normalized to propane (C3). Markers display experimental data while lines display exchange
kinetic model results for each alkane with the same color as the markers. Left: C1–C3
experiment; right: C2–C5 experiment.
Figure 8. An example of paths of alkane hydrogen isotope exchange simulated by the exchange
kinetic model. Temperature = 100
°C.
The trajectories are displayed in solid redlines with the
starting point shown as blue circles. The dotted lines represent equilibrium hydrogen isotope
values.
Figure S1. Torsional Conformers of
n
-butane,
i
-pentane and
n
-pentane.
Figure S2. Equilibrium H isotope fractionation of the H positions, using
n
-pentane (g-g+) and
n
-
butane (g) as examples. Each molecule’s
ε
are referenced to an arbitrary value (6‰ lower than
minimum
β
position of each conformer) for better visualization. Colors are used to separate
hydrogens (with different C–H orientations) on the same carbon.
Figure S3. Equilibrium hydrogen isotope fractionation for alkane pairs vs. temperature. Red line:
Thiagarajan et al. (2020); blue line: this study.
Figure S4. Examples of unequilibrated exchange experiments. All experiments here are
conducted at 100°C.
20