Quantifying the Impact of Parametric
Uncertainty on Automatic Mechanism
Generation for CO
2
Hydrogenation on Ni(111)
Bjarne Kreitz,
∗,†,‡
Khachik Sargsyan,
¶
Emily J. Mazeau,
§
Katr´ın Bl¨ondal,
‡
Richard H. West,
§
Gregor D. Wehinger,
†
Thomas Turek,
†
and C. Franklin
Goldsmith
∗,‡
†Institute of Chemical and Electrochemical Process Engineering, Clausthal University of
Technology, Clausthal-Zellerfeld, Germany
‡School of Engineering, Brown University, Providence, RI, USA
¶Sandia National Laboratories, Livermore, CA 94550, USA
§Department of Chemical Engineering, Northeastern University, Boston, MA 02115, USA
E-mail: kreitz@icvt.tu-clausthal.de(BK); franklin goldsmith@brown.edu(CFG)
Phone: +49 5323 72 2473 (BK); +1 401 863 6468 (CFG)
Abstract
Automatic mechanism generation is used to determine mechanisms for the CO
2
hydrogenation on Ni(111) in a two-stage process, while considering the uncertainty
in energetic parameters systematically. In a coarse stage, all the possible chemistry is
explored with gas-phase products down to the ppb level, while a refined stage discovers
the core methanation submechanism. 5,000 unique mechanisms were generated, which
contain minor perturbations in all parameters. Global uncertainty assessment, global
sensitivity analysis, and degree of rate control analysis are performed to study the effect
1
of this parametric uncertainty on the microkinetic model predictions. Comparison of
the model predictions with experimental data on a Ni/SiO
2
catalyst find a feasible set
of microkinetic mechanisms that are in quantitative agreement with the measured data,
without relying on explicit parameter optimization. Global uncertainty and sensitivity
analyses provide tools to determine the pathways and key factors that control the
methanation activity within the parameter space. Together, these methods reveal that
the degree of rate control approach can be misleading if parametric uncertainty is not
considered. The procedure of considering uncertainties in the automated mechanism
generation is not unique to CO
2
methanation and can be easily extended to other
challenging heterogeneously catalyzed reactions.
Methanation is a promising technology for reducing CO
2
emissions while producing sustain-
able natural gas. From a climate-change perspective, the process is particularly advanta-
geous when excess renewable energy is used to generate the requisite H
2
via water splitting as
part of the Power-to-Gas process.
1,2
However, volatility in renewable energy sources induces
challenges on the transient operation of a catalytic reactor.
3
Given that the net reaction,
CO
2
+ 4 H
2
CH
4
+ 2 H
2
O, is exothermic, ∆H
rxn
(298 K) = −164.7 kJ mol
−1
, transient
operation can lead to undesirable temperature and concentration gradients.
4–6
Accordingly,
an accurate microkinetic mechanism is essential for optimizing reactor performance.
The most commonly used methanation catalyst is Ni, due to its good performance at reason-
able costs.
2,7
Ni(111) has the highest share on a Ni crystal,
8
yet its role in CO
2
methanation
is unresolved, despite extensive research.
9–17
Experiments with Ni/γ-Al
2
O
3
catalysts point to
the higher activity of Ni(111) terrace sites, whereas experiments on Ni/SiO
2
show a higher
activity of Ni(211) steps, which are also considered to be the active site for CO metha-
nation.
18
Lozano et al.
12
combined density functional theory (DFT) calculations using the
BEEF-vdW functional and kinetic Monte Carlo simulations to demonstrate that the Ni(111)
surface is inactive for the CO
2
methanation; instead, they argued that the catalyst converts
the CO
2
to CO in the reverse water-gas shift (RWGS) reaction, CO
2
+ H
2
CO + H
2
O,
2
with the CO
*
2
dissociation being the rate-determining step (RDS). In the study of Vogt et
al.,
10
the four most dominant Ni facets were investigated by DFT with the PBE functional.
They showed in a mean-field microkinetic model that although the Ni(111) facet is not as
active as the open (110) facet, it still exhibited some methanation activity. These authors
identified the dissociation of HCO
*
as the RDS, which is supported by Zhou et al.
15
from
DFT calculations with the same functional. In addition to these models of CO
2
hydrogena-
tion on Ni(111), other studies focused on the (reverse) water-gas shift reaction,
19–21
methanol
synthesis,
22
and formic acid formation.
23
Apart from the general role/activity of the Ni(111)
facet, the dominant reaction network for CO
2
hydrogenation has not been conclusively deter-
mined; the aforementioned computational studies disagree about important intermediates,
pathways, and the RDS.
The microkinetic mechanism for CO
2
hydrogenation on Ni(111) can be developed based
either on surface science experiments assisted by operando methods,
9,10,13,24
or by compu-
tational methods (e.g. DFT).
16,19,25–28
DFT-based microkinetic mechanisms are increas-
ingly common, due to the availability of “black-box” electronic structure codes. Although
modern DFT functionals are reliable for adsorbate thermochemistry and kinetics, these cal-
culations remain computationally expensive. Accordingly, given a computational “budget”,
researchers must prioritize which intermediates and transition states to investigate. This pro-
cess assumes that the researcher knows a priori which intermediates and transition states
will be important. Consequently, the mechanism generation process can be biased by the
developer’s expectations.
29,30
An alternative to DFT-based mechanism development is to use
approximate methods that, while less accurate, are orders of magnitude faster. One such
method is applying linear scaling (LS) relations, which are based on the d-band model.
31
LS
relations can accelerate the procedure,
32
and are often used for catalyst screening.
28,33,34
Al-
though these approximate methods save computational resources, they still require expertise
and intuition to develop the mechanism, and this procedure does not avoid the problem of
incompleteness due to bias. An alternative to constructing mechanisms “by hand” is to use
3
computers to propose and evaluate possible elementary reactions.
35–42
One such approach
is the automatic Reaction Mechanism Generator (RMG) of Green and coworkers.
39,43
Orig-
inally developed for gas-phase pyrolysis, RMG has been expanded to include reactions on
surfaces.
40–42
The omission (intentional or unintentional) of certain reactions can be characterized as mech-
anistic uncertainty. In addition to bias, a second problem with microkinetic mechanism
generation is parametric uncertainty. All of the DFT-derived parameters carry uncertainties
because of the assumptions made in the exchange-correlation functional. For example, the
binding energies are assumed to have an uncertainty of ± 0.3 eV.
44
However, the uncertain-
ties in binding energies for different adsorbates are correlated
45–47
(indeed, some degree of
correlation is implicit in LS, and the BEEF-vdW functional exploits this correlation
48,49
). In
addition to correlation among adsorbate thermochemistry, reaction kinetics are correlated as
well, as exemplified by Brønsted-Evans-Polanyi (BEP) relations.
50
The uncertainty in model
parameters should be propagated to the outputs of the model, e.g. conversion, turnover
frequency (TOF)
44,45,49,51
and to identify the path with the highest occurring frequencies in
a mechanism.
51–53
However, given the large uncertainty in model parameters, some pathways
or intermediates might have been overlooked because of the very complex landscape of the
potential energy surface.
30
Therefore, it is necessary to account for the uncertainty directly
in the mechanism generation procedure to provide an exhaustive analysis of all possible
reactions and intermediates.
The present work aims to combine experiment, theory, and modeling to develop a microki-
netic model for the hydrogenation of CO
2
on Ni(111). Instead of propagating uncertainty
from a final microkinetic model to the simulation outputs, we take a novel approach and
include the uncertainty directly in the mechanism generation procedure in RMG. The au-
tomatic mechanism generation process is repeated 5,000 times, with each new mechanism
including small perturbations in the DFT-derived parameters that can also result in mecha-
4
nisms with different species and reactions. Therefore, we can discover all possible reactions
and intermediates in a vast reaction network. Global sensitivity analysis (GSA) and lo-
cal sensitivity analysis using the degree of rate control (DRC) are used to identify the most
important species and reactions over the whole uncertainty range. The mechanisms are com-
pared against experimental data for a Ni/SiO
2
catalyst in a differential fixed-bed reactor.
This comparison determines a feasible set of microkinetic mechanisms that quantitatively
agree with the experimental data. Combining all of the methods allows us to advance our
understanding of the factors controlling the methanation activity on Ni(111) and to de-
rive a most likely methanation mechanism. Applying the DRC on each unique mechanisms
shows how versatile the DRC can be in a confined uncertainty range and how much more
information on the controlling factors is obtained when global uncertainty is considered.
This study provides an example for the benefit of automated mechanism generation and the
consideration of uncertainties to discover all the possible chemistry.
Materials and Methods
Microkinetic Mechanism Generation
Automated mechanism generation was performed with RMG (version 3.0).
43,54
A detailed
explanation of the RMG software can be found in the work of Gao et al.
39
and the exten-
sion to heterogeneously catalyzed reactions in the publication from Goldsmith and West.
40
Only a brief explanation of the key features important for this work is provided. Elemen-
tary reactions are grouped according to reaction families, which are templates that convert
the chemical graphs of reactants into products. For each proposed species and elementary
reaction, RMG must provide thermodynamic properties and a rate coefficient, respectively.
RMG combines a database of precompiled values, but it can supplement these databases
with rules for predicting the properties for novel species/reactions. The total species popu-
lation is divided into two groups: core and edge. Core species are essential to the mechanism;
5
