Ab initio coverage-dependent microkinetic modeling of benzene hydrogenation on Pd(111)

Maarten K. Sabbe *a, Gonzalo Canduela-Rodriguez ab, Jean-François Joly b, Marie-Françoise Reyniers a and Guy B. Marin a
aLaboratory for Chemical Technology, Universiteit Gent, Technologiepark 914, 9052 Gent, Belgium. E-mail: maarten.sabbe@ugent.be
bIFP Energies nouvelles, BP 3 - 69360 Solaize, France

Received 15th May 2017 , Accepted 11th September 2017

First published on 12th September 2017


The effect of hydrogen coverage on the kinetics of benzene hydrogenation on Pd(111) has been investigated with optPBE-vdW density functional theory calculations and a coverage-dependent microkinetic model. The dominant reaction path consists of the consecutive hydrogenation of carbon atoms located in ortho positions relative to the previously hydrogenated carbon atom, independent of the hydrogen coverage. Increasing the hydrogen coverage destabilizes all surface species, which leads to weaker adsorption and increased rate coefficients for the hydrogenation steps due to stronger destabilization of reactants than transition states. The catalytic activities simulated using the constructed coverage-dependent microkinetic model exceed those obtained using a low-coverage microkinetic model by several orders of magnitude and are comparable to experimentally observed activities. The rate coefficients to which the global rate is most sensitive depend on the reaction conditions and differ from those calculated using low coverage kinetics. Therefore, properly accounting for coverage dependence on the kinetics and thermodynamics of catalytic hydrogenation reactions is not only required for an accurate DFT-based prediction of the catalytic activity but also for a correct understanding of the reaction mechanism.


1. Introduction

Density functional theory (DFT) based microkinetic modelling has evolved rapidly in the past decades. The combination of DFT and microkinetic modelling, in which catalytic reactions are described in terms of elementary reaction steps,1 allows the exploration and unravelling of reaction mechanisms. Furthermore, the use of DFT-based kinetics in microkinetic models avoids the error compensation that can occur during the standard approach of estimating the parameters in the microkinetic model. The estimation of the kinetic and thermodynamic parameters of a microkinetic model can obscure deficiencies in the microkinetic model, possibly leading to an incorrect chemical understanding of the process. However, despite the rapid advances in computational chemistry in the past decades and the successes in structure determination and qualitative assessment of reactivity, the combination of DFT and microkinetics rarely succeeds to quantitatively predict the activity of heterogeneous catalysts. There are a myriad of possible reasons and unaccounted imperfections that can be at the root of this,2 but this paper focuses on the necessity of performing the DFT calculations at a relevant coverage to reliably predict the activity of metal-catalyzed processes, with benzene hydrogenation over Pd(111) as the reaction under study.

For metal-catalyzed heterogeneous catalysis, density functional theory calculations have for the major part been oriented towards describing the metal–adsorbate bond in its purest form, avoiding lateral interactions between neighbouring adsorbates. Using periodic boundary conditions, as is standard practice for calculations on low Miller-index surfaces, this means that the lateral interactions are minimized using the largest possible supercell that is computationally affordable. These catalyst models will typically represent low to medium coverages, which is well-suited for many low-pressure processes. However, association reactions that are performed at higher pressures, such as hydrogenation reactions, are typically associated with higher surface coverages. Since adsorbate–adsorbate interactions can affect chemical kinetics to a large extent,3–8 kinetic models based on DFT calculated kinetics at low coverages are often not able to predict experimental turnover frequencies or the experimentally observed kinetic trends.

To account for the effect of reaction conditions on the strength of adsorbate–adsorbate interactions and its effect on the reaction kinetics means that coverage dependence of the kinetic and thermodynamic parameters needs to be included. To this end, semi-empirical approaches have already been developed, many to explicitly account for neighboring adsorbates in kinetic Monte Carlo methods. These methods are typically based on the bond order conservation/unity bond index quadratic exponential potential method (UBI-QEP)9 and have been implemented in the studies of, e.g., the groups of Neurock,6 Vlachos,3,10 Reuter,11 Mavrikakis12 and Norskov.7 These models typically combine theoretical and experimental results,10–12 including regression of parameters to experimental data.3 Such models are limited to reactions between small compounds, such as H2 oxidation and ethylene hydrogenation, since the application of UBI-QEP methods to predict binding of larger adsorbates with catalyst surfaces is not straightforward. Also, more complicated cluster expansion approaches, which expand adsorbate energies on a general basis of pair, triplet, and higher order adsorbate–adsorbate interactions, have only been applied yet to adsorption/desorption processes or catalytic reactions of very small molecules as NO.13 For catalytic reactions of larger reactants, kinetic Monte Carlo methods prove already highly complicated, and accounting for lateral interactions in the kMC simulation of the reaction of larger reactants, with the associated larger reaction mechanisms, is virtually impossible. Therefore, accurate modelling of coverage dependence for reactions of larger compounds can better be performed using DFT methods directly, and the results incorporated in standard microkinetic models. In order to calculate coverage dependence of the thermodynamic and kinetic parameters with DFT methods, the first step is to evaluate the range of surface coverages that can be expected on the catalyst surface under relevant reaction conditions based on calculation of surface Gibbs free energies.2,14,15 In a second step, the reaction kinetics and thermodynamics can be calculated for the expected range of surface coverages under reaction conditions. Ideally, all possible coverages should be evaluated, but this is computationally not feasible. A range of coverages can, however, be evaluated and applied to construct a coverage-dependent microkinetic model, by calculating activation energies and entropies as a function of the fractional coverage. This explicit coverage-dependent model can be solved to predict production rates and surface concentrations under reaction conditions.

In this work, the target reaction for coverage-dependent modelling is catalytic benzene hydrogenation over Pd(111). For this reaction, there is still discussion about some aspects of the reaction mechanism, but at the same time, there are sufficient experimental data available to compare the predicted reaction kinetics to. The available experimental and computational research into benzene hydrogenation over Pd is summarized in the background section. To the best of our knowledge, no such coverage-dependent model has been reported yet for molecules the size of benzene that is based on ab initio results only.

This work therefore aims at constructing a coverage-dependent microkinetic model for benzene hydrogenation on Pd(111) based on periodic DFT calculations, in which explicit coverage effects on both adsorption and hydrogenation kinetics are included. To this end, the dominant hydrogenation path is identified and evaluated at various, increasing hydrogen coverages on the surface. The results are used to construct a coverage-dependent microkinetic model that is integrated to predict catalytic activities. The microkinetic results are compared to experimental results from the literature.

Benzene hydrogenation: background

Hydrogenation of aromatic molecules is important in the refining industry to improve the quality of diesel fuels and to produce clean fuels that meet the stringent legislation on the concentration of carcinogenic compounds in fuels.16 The most straightforward way to model these processes is to study the hydrogenation of benzene, the least complex aromatic reactant. Benzene hydrogenation is also important to the petrochemical industry as the first step in the synthesis of caprolactam and adipic acid, which are monomers for the production of nylon 6 and 6-6 polymers.17 A commonly used hydrogenation catalyst is Pd, yet its high cost is a strong driving force to optimize hydrogenation catalysts further to higher activities/selectivities at a lower cost. Alloying Pd with other metals is an evident method to reduce the cost, with opportunities to improve the catalyst properties at the same time. However, to enable the development of Pd-based bimetallic hydrogenation catalysts based on fundamental understanding, the mechanism on Pd needs to be well understood first. However, despite the large number of experimental and theoretical studies on benzene hydrogenation, there is no consensus on the mechanism of this complex reaction network.8 Microkinetic modelling, in which catalytic reactions are described in terms of elementary surface steps, can help to explain the mechanism.1 Most details about these elementary steps can be obtained using computational chemistry techniques such as density functional theory (DFT), which can be very helpful to unravel the complexities of catalytic reaction mechanisms.

The general agreement in the literature is that the catalytic benzene hydrogenation proceeds via the Langmuir–Hinshelwood–Hougen–Watson mechanism involving reaction between chemisorbed hydrogen atoms and chemisorbed hydrocarbon species. These findings are based on experimental work on the main hydrogenation catalysts Ni,18–25 Pd,24,26–28 and Pt,24,28–35 with only a few authors postulating that the Eley–Rideal mechanism is followed, with gas-phase hydrogen reacting with chemisorbed benzene on Ni and Cu.36,37 Despite this consensus, a wide range of kinetic models has been proposed.

Competition between benzene and hydrogen for the active sites has been both claimed24,25,28 and denied.18,19,22,23,38 As for the rate determining step, the first hydrogen addition is often considered as the rate determining step on Ni,20,21,38 Pd,38 and Pt.30,34 On Pt, also the fourth39 and fifth hydrogen addition40 have been proposed as rate determining steps, and on Ni the assumption has been made that all six sequential hydrogenation steps have the same rate coefficient.18,21 Some models incorporate the dehydrogenation of benzene, which can explain the formation of carbonaceous species30 and the subsequent deactivation caused by a decrease in active sites.27 However, benzene dehydrogenation is neither thermodynamically nor kinetically favorable compared to hydrogenation.41 The presence of adsorbed cyclohexene during benzene hydrogenation experiments has been reported, particularly on Pt,31,32,42–44 although the formation of gas-phase cyclohexene is thermodynamically less favorable than the formation of the major product cyclohexane.42,45 Inhibition of the reaction by adsorbed cyclohexane can be discarded46 because its surface concentration is usually assumed to be negligible, since cyclohexane is typically expected to desorb fast and irreversibly.47 The catalytic activity reaches a maximum as a function of temperature at about ∼500 K for Pd catalysts,26,28 explained by the balance between increasing rate coefficients and decreasing reactant surface concentrations with temperature. The partial reaction order with respect to hydrogen is usually in the range from 0.5 to 3,26 while it is close to zero for benzene.24,26 This could indicate that the fractional surface coverage of benzene is higher than that of hydrogen.

There are few theoretical studies into benzene hydrogenation,8,47–49 and these usually consider low to moderate reactant coverages. Based on DFT calculations, Mittendorfer and Hafner48 showed that the electronic reaction barriers for benzene hydrogenation on Ni(111) via the Langmuir–Hinshelwood mechanism are much lower than for the Eley–Rideal mechanism. On Pd(111), Morin et al.49 derived from thermodynamic considerations that the reaction path follows the Horiuti–Polanyi mechanism50 with successive hydrogenation of neighboring carbon atoms, i.e. in ortho positions relative to the previously hydrogenated carbon atom (ortho path). Based on the Evans–Polanyi relationship, the first hydrogen addition was identified as the rate limiting step. The same step was proposed on Ni(111) from activation barriers, in line with what was previously proposed in the literature, because the aromaticity of benzene is not broken upon adsorption.20,21,38 On Pt(111), initially the reaction path involving the hydrogenation of carbon atoms in the meta position (for the first three hydrogenation steps) was found to be preferred on Pt22 clusters,47,49 with the fifth reaction having the largest activation energy.47 However, the ortho path was later identified based on regression to experimental data51 and confirmed by more recent DFT work.8 The ortho path also agrees with the observations of cyclohexene formation during both benzene hydrogenation35 and cyclohexane dehydrogenation4 on Pt catalysts.

2. Methodology

2.1. Electronic energy calculations

The Vienna Ab initio Simulation Package (VASP 5.3.3)52–55 is used to perform periodic density functional theory (DFT) calculations. The projector augmented wave method (PAW)56,57 is used to describe the valence electrons, with a 400 eV cut-off energy for the plane waves. The first-order Methfessel Paxton method58 is used to describe partial electronic occupancies close to the Fermi level with a smearing width of 0.3 eV and Brillouin-zone integration is performed on a 5 × 5 × 1 Monkhorst–Pack grid.59 To describe the electron–electron interactions, two exchange–correlation functionals are used: for the pre-screening of the full reaction network the generalized gradient Perdew–Wang PW91 is used,60,61 while more accurate calculations on the dominant path are performed using the optPBE-vdW functional developed by Rydberg and Dion62,63 and implemented in VASP by Klimeš et al.64,65 The latter functional has been shown to provide adsorption energies and geometries for benzene adsorption on transition metal surfaces in line with experimental observations, thanks particularly to an improved description of the nonlocal nature of the electron correlation and in particular van der Waals interactions.66,67 Using these settings, the calculated magnetization of the bulk is below 10−4μB and the adsorption energy of benzene is within 6 kJ mol−1 of a spin polarized calculation. This agreement is judged sufficient to justify spin unpolarized calculations in this work.

The catalyst is modeled as a Pd(111) surface, representing the most abundant surface facet on larger catalyst particles. Pd hydride or subsurface hydrogen has not been considered since earlier research indicates this not to be stable under typical industrial hydrogenation conditions of 450 K and hydrogen pressures up to hundreds of bar.14,15 The computational model is a 4-layered slab with a (3 × 3) unit cell, which is sufficiently accurate to describe surface reactions (see Fig. 1).49 The associated fractional coverage, with 1 benzene per 9 Pd surface atoms, represents a moderate to high benzene coverage as will be found under benzene hydrogenation conditions according to our previously constructed phase diagrams.15 The two bottom Pd layers are fixed to the optimized bulk structure (calculated lattice parameter a = 395.7 pm), while the two top Pd layers and the adsorbates on the surface are fully relaxed during the geometry optimization. The unit cell is repeated to model an infinite surface, and periodic boundary conditions are applied. An 11 Å vacuum gap is present between the consecutive slabs, containing a dipole layer in order to avoid dipole interactions between periodic images. Very strict convergence criteria are applied in order to avoid imaginary frequencies in the subsequent frequency analyses, amounting to 10−8 eV on the energy in the electronic minimization and 0.015 eV Å−1 for the maximum force in the geometry optimization. Vibrational analyses, which are required to calculate thermodynamic properties, are performed in the harmonic oscillator approximation. VASP numerically evaluates the Hessian matrix with a step size of 0.015 Å. The vibrational frequencies are calculated for all atoms that are relaxed in the geometry optimization, i.e., all atoms of the two top Pd layers of the slab and the adsorbates.


image file: c7cy00962c-f1.tif
Fig. 1 Four-layered 3 × 3 supercells used in the periodic DFT calculations. Top view (left, only 2 layers shown) and side view (right) (Pt: dark grey; C: light grey; H: white).

In previous work15 we have shown that under typical hydrogenation conditions the surface of Pd(111) is populated with higher hydrogen coverages than usually applied in DFT studies: up to θH = 0.89 under typical hydrogenation conditions. Therefore, in this work the effect of hydrogen coverage on the thermodynamics and kinetics of benzene hydrogenation is evaluated. To this end, adsorption and reaction steps are studied at four different hydrogen coverages (θH,initial). The hydrogen atoms are coadsorbed with the hydrocarbon intermediates along the reaction path. The hydrocarbon coverage itself is fixed to a single molecule per (3 × 3) unit cell. With the coverage defined as the number of adsorbed molecules per surface Pd atom, this corresponds to a coverage θBHi of 0.11. At a coverage of 1.64 × 1018 molecules per m2, this is 70% of the monolayer coverage of 2.3 × 1018 molecules per m2 experimentally reported on Pt(111) at 300 K.5 Our earlier evaluations have shown that under typical hydrogenation conditions, the studied fractional benzene coverage of θBHi = 0.11 is the most probable coverage on Pd(111).15

The four evaluated hydrogen coverages pertain to θH,initial = 0, 0.11, 0.44 and 0.67, with the subscript ‘initial’ pointing at the hydrogen coverage coadsorbed with the reactant, before one of the hydrogen atoms adds to the reactant. For the first case of θH,initial = 0, it is considered that the reacting hydrocarbon and hydrogen atom are each adsorbed in different unit cells (θBHi = 0.11 and θH = 0.11 in separate cells without coadsorption). The other three coverages θH,initial = 0.11, 0.44 and 0.66 consider 1, 4, and 6 hydrogen atoms coadsorbed in the (3 × 3) unit cell with a single hydrocarbon (θBHi = 0.11). The hydrocarbon is then hydrogenated by one of these coadsorbed hydrogen atoms. The coadsorbed hydrogen atoms are assumed to always adsorb on the sites that correspond to the most stable adsorption sites for coadsorption with benzene, regardless of the degree of saturation of the hydrogenated intermediate. These sites correspond to those calculated in a previous study for coadsorption of benzene with different hydrogen coverages.15

Transition states are obtained using the following approach. First, the minimum energy path between the optimized reactant and product state is calculated with the nudged elastic band method.68 The preliminary saddle point geometry resulting from this NEB is further used as the initial guess for the dimer method.69 At a high hydrogen coverage, the optimization of transition states can be difficult to converge due to the increased number of atoms and the steric interactions between the coadsorbed species. In order to accelerate convergence, the transition state found for the lowest hydrogen coverage θH = 0 is used as the initial geometry for the optimization at increased hydrogen coverages, after adding the proper amount of hydrogen atoms in the unit cell. This approach yields the same transition state as optimization from scratch, starting with an NEB calculation, as has been checked for the first hydrogenation step at θH,initial = 0.44.

2.2. Thermodynamic and kinetic parameters

Standard statistical thermodynamic expressions, in the harmonic oscillator approximation, and conventional transition state theory are used to calculate kinetic and thermodynamic parameters as a function of temperature, as described in detail in section S1 of the ESI. In the following section, the salient aspects related to the modeling of coverage dependence are discussed.

For the calculation of pre-exponential factors, all chemisorbed intermediates are assumed to be immobile under reactive conditions; even for hydrogen, typically one of the most mobile species, the electronic barrier for hydrogen diffusion between the most stable adsorption sites is above 15 kJ mol−1 at half the saturation coverage. Therefore, hydrogen is, even at the hydrogenation temperature of 450 K (RT = 3.7 kJ mol−1), below the threshold for free mobility. Hydrogen is certainly mobile in an activated sense, but does not correspond to free translation, and its exact partition function will rather approach the harmonic oscillator approximation than free translation. Therefore, all chemisorbed surface species are assumed to be immobile for the calculation of the partition function, and only vibrational contributions are considered for surface species, including the translational and rotational motions of surface species. For physisorbed cyclohexane, the two imaginary frequencies that could not be avoided were replaced by frequencies of 25 cm−1, a typical value for this kind of large-amplitude motion. For gas-phase species, free rotation and translation is considered next to the vibrational contributions.

The adsorption rate coefficient is described as the product of the incident Hertz–Knudsen molecular flux F and the sticking probability on a clear surface s0, as shown in eqn (1). Adsorption is considered as a non-activated process for all species, with the sticking probability s0 assumed to be 1 for all adsorbates. Desorption rate coefficients are obtained using the thermodynamic equilibrium coefficient from DFT adsorption enthalpies and entropies, as shown in eqn (2).

 
image file: c7cy00962c-t1.tif(1)
 
image file: c7cy00962c-t2.tif(2)

In eqn (1), F is the incident molar flux (s−1), p the pressure in Pa, Nt is the number of active sites per m2 (1.64 × 1018 m−2), m is the molecular mass (kg), kB is the Boltzmann constant and T is the temperature.

Surface rate coefficients (k in s−1) are described with an Arrhenius expression, using the activation energies (Ea) and pre-exponential factors (A) calculated as explained in ESI section S1. To ensure thermodynamic consistency in the microkinetic model,1 the kinetic parameters for the reverse reaction are calculated from those for the forward reaction:

 
Ea,rev = Ea,for − ΔrH0(3)
 
image file: c7cy00962c-t3.tif(4)
where Ea,for and Ea,rev are the forward and reverse activation energies, Afor and Arev the forward and reverse pre-exponential factors, and ΔrH0 and ΔrS0 are the standard reaction enthalpy and entropy. The procedure for calculating these parameters is explained in detail in section S1 of the ESI. The consistency of the calculated equilibrium coefficients has been verified: the product of the equilibrium coefficients of all reactions along the dominant path is exactly the same as the calculated equilibrium coefficient of the gas phase reaction between benzene and cyclohexane.

For the calculation of kinetic and thermodynamic parameters at higher hydrogen coverage, the hydrogen coverage coadsorbed with the reactants (θH,initial) is assumed to remain constant along the path. After every hydrogenation step, in which one of the coadsorbed hydrogen atoms adds to the hydrocarbon, the hydrogen coverage coadsorbed with the product (θH,final) will be lower than θH,initial (see Fig. 2 and expression (5)).

 
image file: c7cy00962c-t4.tif(5)
 
image file: c7cy00962c-t5.tif(6)


image file: c7cy00962c-f2.tif
Fig. 2 Schematic representation of a hydrogenation step with subsequent hydrogen addition to restore the initial hydrogen coverage θH,initial. The coverage of hydrogen coadsorbed with the reactants (θH,initial) is higher than with products (θH,final) because one hydrogen has been added to the hydrocarbon. Therefore, one hydrogen atom is adsorbed on the surface covered with θH,final and the product hydrocarbon to restore the initial θH,initial coverage.

Therefore, the hydrogen coverage in the product unit cell is restored to θH,initial by adsorbing an additional hydrogen atom in the unit cell (see eqn (6) and Fig. 2). This yields again the initial coverage θH,initial for the next hydrogenation step. This approach leads to six separate hydrogen adsorption steps, after every hydrogenation step. The hydrogen adsorption enthalpy and entropy for every hydrogen adsorption step are calculated from eqn (6). Implementation of the kinetics in a mean-field microkinetic model, however, requires a unique value for the hydrogen adsorption enthalpy, and not six different values (one for every hydrogenation step). Therefore, the average value of the hydrogen adsorption enthalpy is used in the kinetic model.

For an accurate simulation of activities, the accuracy of the overall equilibrium from benzene to cyclohexane, as well as to the side products cyclohexene and cyclohexadiene, is of prime importance. However, the gas-phase equilibrium coefficient as calculated by DFT does not match the experiment-based equilibrium coefficient that can be calculated from data in the National Institute of Standards and Technology database (NIST).70 This is mainly caused by an overestimation in the reaction enthalpy: for benzene to cyclohexane the DFT reaction enthalpy overestimates the experimental value by 22 and 41 kJ mol−1 for the PW91 and optPBE-vdW functional, respectively. The difference in reaction entropy is much smaller, below 3 J mol−1 K−1 for benzene to cyclohexane. To correct for the deviation of the gas phase equilibrium from experiment, the adsorption enthalpies and entropies for the three product species, i.e., 1,3-cyclohexadiene, cyclohexene and cyclohexane, are adapted in order to match the thermodynamic gas phase equilibrium as calculated from the NIST database.8 The adsorption and desorption enthalpies and enthalpies for the three product species shown in this paper are the adapted values unless noted differently; the original ab initio values are given in the ESI as well. Details of the calculation can be found in ESI section S2.

To enable the application of the kinetic and thermodynamic data calculated at four different hydrogen coverages θH,initial into a microkinetic model, a continuous function of these parameters as a function of coverage is required. This allows the evaluation of the kinetic parameters at the effective coverage resulting from the microkinetic simulation. In order to include an explicit function of the thermodynamic and kinetic parameters as a function of coverage in the microkinetic simulation, a linear regression of the thermodynamic and kinetic parameters as a function of the total coverage θtotal has been made (see eqn (7)). The total fractional coverage, accounting for the coverage of all surface species, is taken as θtotal = 1 − θ*, i.e., 1 minus the fraction of free sites.

 
L = total + b(7)

In this equation, L is a parameter in the Arrhenius equation or adsorption equilibrium coefficient (Ea, ΔS, ΔHads, ΔSads). The constructed coverage-dependent functions can then be implemented into the microkinetic code to provide direct coupling between the simulated coverage and the coverage at which the kinetics and thermodynamics are evaluated. The coverage-dependent rate coefficients are evaluated at the total coverage that results from the microkinetic simulation, which makes solving a coverage-dependent microkinetic model an iterative procedure.

Microkinetic model and reactor simulation. Catalytic activities are simulated, solving the microkinetic model as a function of reaction conditions in a reactor model. The reactor model allows a realistic build-up of product pressures, without which the predicted turnover frequency would be too high since no readsorption or equilibrium setting is allowed for in the absence of product pressures.

The microkinetic rate equations are as follows. The adsorption rate rj of species i is calculated according to:

 
image file: c7cy00962c-t6.tif(8)
with the fractional coverages θ squared for dissociative adsorption of H2. The rate rj of every surface reaction, eqn (9), is calculated using eqn (10).
 
image file: c7cy00962c-t7.tif(9)
 
image file: c7cy00962c-t8.tif(10)

The sum of all fractional surface coverages amounts to one:

 
image file: c7cy00962c-t9.tif(11)

There are no two separate site balances used for hydrogen and the hydrocarbons, which is sometimes done to model hydrogenation processes, because this will prohibit the simulation of the negative partial reaction orders in benzene that are observed experimentally at lower temperatures.1,26

The net production rate Ri of surface species i (molecules per site per s) should be 0 in the pseudo-stationary state approximation:

 
image file: c7cy00962c-t10.tif(12)

The pseudo-steady-state surface concentrations are obtained transiently by solving the differential equations of eqn (12) till a steady state is reached, using the LSODA integration routines from the ODEPACK Fortran library as implemented in the Python package SciPy.

The resulting microkinetic model is embedded in a reactor model that algebraically solves the outlet flow rates of component i in the mass balance for a differentially operated reactor:71

 
Fi,inFi,out + WCtRi = 0(13)

In this equation, Fi is the molar flow of component i (mol s−1), W is the catalyst mass (1.28 × 10−3 kg), Ct is the concentration of active sites (9.6 × 10−2 molactivesite kgcat−1, as obtained from CO chemisorption on Pd powder26), and Ri is the net formation rate of species I expressed as a turnover frequency (s−1). The conversion is in all simulated cases below 5%, which justifies the assumption of a differentially operated reactor. At these low conversions, the equations can also be used to simulate the tubular reactors on which some of the experimental reference data were determined.

The partial reaction order is evaluated by varying its inlet pressure and keeping all other inlet pressures constant apart from the pressure of the inert (N2). The sensitivity of the catalytic activity to kinetic and thermodynamic parameters of the elementary reaction steps is quantified by means of the sensitivity coefficient. To this end, the pre-exponential factor of all reactions are varied one by one, while the kinetic parameters of all other reactions are kept constant; see eqn (14) below. The sensitivity coefficient XSA,j expresses the relative change in the turnover frequency in proportion to the relative change in the pre-exponential factor,

 
image file: c7cy00962c-t11.tif(14)
with TOF as the turnover frequency calculated with the original parameters, TOF′ the turnover frequency calculated with the modified pre-exponential factor of step j, Aj the original pre-exponential factor and image file: c7cy00962c-t12.tif the modified pre-exponential factor, taken as 1.01 times Ai. The thermodynamics are not changed during the sensitivity analysis, meaning that the forward and the reverse pre-exponential factors of a reaction are always changed simultaneously.

3. Results and discussion

Periodic DFT calculations have been performed to study the catalytic hydrogenation of benzene to cyclohexane on Pd(111). First, the results at low coverages are presented, i.e. for the full reaction network without coadsorbed hydrogen in the DFT calculations (“low coverage”), and the dominant path is identified. Next, the effect of hydrogen coverage on the adsorption and hydrogenation reactions is studied on this dominant path. The DFT-calculated kinetics and thermodynamics are further used to construct a coverage-dependent microkinetic model that predicts catalytic activities under actual reaction conditions.

3.1. Kinetic analysis at low coverage

There are multiple adsorption sites on the Pd(111) surface on which the reactants benzene and hydrogen can adsorb. In principle, the reactions between benzene and hydrogen adsorbed at all possible combinations of adsorption sites would have to be considered. However, to limit the number of DFT calculations, benzene is assumed to adsorb only at the two most stable adsorption sites, bridge and hollow-hcp, as determined earlier.15 Furthermore, it is assumed that the partially hydrogenated hydrocarbon intermediates remain adsorbed on the same site as benzene, and no diffusion to another adsorption site occurs, resulting in two distinct kinetic models, a ‘bridge’ and a ‘hollow’ model. Even with this approximation, there are, for every benzene adsorption site, twenty different hydrogenation steps from benzene to cyclohexane via 13 possible reactions paths (see Fig. 3).
image file: c7cy00962c-f3.tif
Fig. 3 Full reaction network of the catalytic hydrogenation of benzene to cyclohexane, with numbering of the hydrogenation steps and with the identified dominant path indicated with green arrows. The following nomenclature is used for the surface species: B: benzene, BH for monohydrobenzene, 13CHD for 1,3-cyclohexadiene, 13DHB for 1,3-dihydrobenzene, 14CHD: 1,4-cyclohexadiene, 123THB: 1,2,3-trihydrobenzene, 124THB: 1,2,4-trihydrobenzene, 135THB: 1,3,5-trihydrobenzene, CHE: cyclohexene, 1235THB: 1,2,3,5-tetrahydrobenzene, 1245THB: 1,2,4,5-tetrahydrobenzene, C-hexyl: cyclohexyl, CHA: cyclohexane.

The aim is to firstly identify the dominant path in the full reaction network at low coverages using the PW91 functional. In this way, the calculations are sufficiently straightforward to enable the calculation of the kinetics and thermodynamics for the whole reaction network on both the bridge and the hollow sites. The low surface coverage involves the reacting hydrocarbon and hydrogen to be adsorbed in different unit cells for the DFT calculation of the thermodynamic and kinetic parameters, which avoids repulsive interactions between hydrocarbon and hydrogen species. To determine the dominant path in the reaction network, gas-phase cyclohexane is assumed to be the only reaction product, which simplifies the microkinetic calculations. As explained in the methodology section, the cyclohexane adsorption equilibrium is adapted in order to set the gas phase equilibrium between benzene and cyclohexane equal to the experimentally-based equilibrium coefficient derived from NIST data.70

The adsorption of a single hydrogen atom in the applied (3 × 3) unit cell corresponds to θH = 0.11, 1 atom per 9 surface Pd atoms, which is 1/9th of an experimentally observed monolayer coverage, corresponding to 1 H atom per surface Pd atom (saturation coverage of 1.47 × 1019 atoms per m2 (ref. 72)). At θH = 0.11, the calculated adsorption enthalpy of hydrogen amounts to −125 kJ molH2−1 (at 450 K) for the most stable hollow fcc site, shown in the right panel of Fig. 4. This adsorption energy is more negative than the calculated value in previous theoretical studies at lower coverages14,73 and also 40 kJ mol−1 more negative than the experimental value for coverages up to θH = 0.5.74


image file: c7cy00962c-f4.tif
Fig. 4 Preferred adsorption sites for benzene, (a) bridge(30) and (b) hollow-hcp(0), and for hydrogen on the (c) hollow fcc surface site of Pd(111) (shown in red for easier identification). The second Pd layer is shown in grey for clarity.

For benzene, the two most stable adsorption sites are the 4-fold bridge(30) and 3-fold hollow-hcp(0) (see Fig. 4a and b), in agreement with previous studies on Pd(111)15,75 and other (111) transition metals.76–78 The calculated adsorption enthalpies on these sites are −113 and −96 kJ mol−1, respectively, at 450 K in a (3 × 3) unit cell. The associated fractional coverage, with 1 benzene per 9 Pd surface atoms, formally written as θB = 0.11, represents a moderate to high benzene coverage, with Tysoe et al. reporting a saturation coverage at θB = 0.13 for planar adsorbed benzene and θB = 0.22 for tilted benzene.79 The adsorption of both bridge and hollow adsorbates is weaker than the experimental range of −140 to −196 kJ mol−1,80 which may indicate that the gradient-corrected PW91 functional is not well-suited to properly describe benzene adsorption. We have also shown in a previous study that the hollow adsorbate is more stable than the bridge adsorbate when coadsorbed with hydrogen.15 Therefore, the reaction network is not only studied starting from bridge-adsorbed benzene, but also starting from hollow-adsorbed benzene, resulting in two reaction networks defined as the ‘bridge’ and ‘hollow’ networks.

The calculated kinetic and thermodynamic parameters of every elementary reaction step in both bridge and hollow networks are given in Tables S4 and S5 in the ESI. The bridge and hollow microkinetic models based on these data are solved under actual reaction conditions (see section 0) using the conditions of Chou and Vannice (W/FB,in = 62.5 kgcat s mol−1; FB,in = 2.1 × 10−5 mol s−1; FH2,in = 2.8 × 10−4 mol s−1, pB,in = 6.7 × 10−2 bar; pH2,in = 0.9 bar).26 The simulated turnover frequency (TOF), defined as molecules of benzene consumed per active site per second (s−1), is larger for the hollow network than for the bridge network (up to 30 times larger at 400 K and 1.5 larger at 550 K) (see Fig. 5a). The TOF for both networks follow the same trend with temperature with a maximum at 550 K, while experimentally this is usually observed at 500 K.26,28 This maximum is clearly linked with the decrease in benzene and hydrogen coverage with temperature (see Fig. 5b). Below 550 K, the decrease in reactant coverage, which should decrease the TOF, is compensated by the increase of the rate coefficients with temperature. Above 550 K the coverage decrease is no longer compensated by the increase in rate coefficient with temperature.39 Because of the higher reactivity of the benzene “hollow” network, only this network will be further discussed below.


image file: c7cy00962c-f5.tif
Fig. 5 (a) Simulated turnover frequency (TOF) (s−1) as a function of temperature (K) for the bridge network (full line) and the hollow network (dotted line). (b) Fractional coverages of hydrogen (black), benzene (grey) and free sites (red) as a function of temperature, obtained for the simulation of the bridge network (full lines) and hollow network (dotted lines). The full and dotted red lines overlap (W/FB,in = 62.5 kgcat s mol−1; FB,in = 2.1 × 10−5 mol s−1; FH2,in = 2.8 × 10−4 mol s−1, pB,in = 6.7 × 10−2 bar; pH2,in = 0.9 bar). Simulations are performed using the kinetic parameters in Table S4 for the network at the hollow site and those in Table S5 at the bridge site.

The calculated turnover frequency at 413 K (TOF = 1.4 × 10−6 s−1) is more than three orders of magnitude lower than the experimental value of 6 × 10−3 s−1 obtained by Chou and Vannice on Pd powder under the same conditions.26 The calculated partial reaction orders for benzene and hydrogen at 450 K are ∼0.2 and ∼4.2, respectively. The value for benzene is close to the experimentally observed zero order;26,81,82 however, the order in hydrogen is larger than the experimental observations, which range between 0.5 and 3.5.26

Based on the net rates resulting from the microkinetic simulations, the dominant reaction path can be identified. The dominant path, shown in Fig. 3, proceeds from benzene (B) via monohydrobenzene (BH), 1,3-cyclohexadiene (13CHD), 1,2,3-trihydrobenzene (123THB), cyclohexene (CHE) and cyclohexyl (C-hexyl) to cyclohexane (CHA). This path corresponds to consecutive additions of hydrogen atoms to the carbon atom in the ortho position relative to the previously hydrogenated carbon, and is therefore called the “ortho path”. This dominant ortho path has also been suggested on Pd(111) by Morin et al.49 on a thermodynamic basis, and by Tetenyi and Paal42 from the hydrogenation of 14C-labelled benzene and cyclohexadiene. The dominant path accounts for over 90% of the overall benzene conversion over the entire temperature range studied (see Fig. S1 in the ESI). Therefore, in the remainder of this work, only the dominant ortho path of the “hollow” network will be considered.

Along this dominant ortho path, the steps to which the overall TOF is most sensitive are the third H addition step for T < 500 K and the fifth H addition step for T > 500 K. At 450 K, the sensitivity coefficient for the third hydrogen addition (reaction d1 in Fig. 3) is XSA,d1 = 0.8, and for the fifth reaction (reaction f1 in Fig. 3), XSA,f1 = 0.2 (see Table S6 of the ESI for more details). At 550 K, the sensitivity coefficient for the fifth H addition increases to a value of XSA,f1 = 0.8 at 550 K and becomes the highest sensitivity coefficient. All the other sensitivity coefficients along the dominant path are zero. This in contrast to what one would expect based on the activation barriers, with the highest barrier along the dominant path present for the first hydrogenation step B* + H* → BH* + * (reaction b in Fig. 3). The first H addition being rate determining has also been suggested by other authors, based on reaction equilibrium coefficients derived from DFT calculations,49 or based on the theoretical consideration that because the aromaticity of benzene is not broken on Pd(111), it should break during the first H addition step.20,21,38 However, the sensitivity coefficient for the first hydrogenation step calculated in this work is zero. The reason why the overall rate is not sensitive to the first reaction is due to the fact that the fractional coverage of benzene, the reactant of the first H addition, is larger than those of the other hydrocarbon intermediates on the surface (see Fig. 5b), which results in an actual forward rate of the first hydrogenation step that is larger than those of the other steps.

The large underestimation of the experimentally observed turnover frequency by the simulation could be caused by insufficient benzene coverage, which strongly decreases with temperature due to too weak adsorption (see Fig. 5b). Every attempt to obtain a quantitative agreement between simulated and experimental turnover frequency values should start with a correct prediction of the adsorption equilibria. A functional that better accounts for non-local effects, such as the optPBE-vdW functional, is expected to describe benzene adsorption better than the PW91 functional used to screen the full bridge and hollow networks. Therefore, the hollow dominant path has been evaluated at low coverages using the optPBE-vdW functional. The dominant path obtained with the optPBE-vdW functional remains the same as the one calculated with the PW91 functional; every reaction of the dominant path has an activation energy that is at least 15 kJ mol−1 lower than the values calculated for the alternative reactions (see Table S7 in the ESI). A benzene adsorption enthalpy ΔHads of −163 kJ mol−1 is calculated at 450 K and θB = 0.11 with the optPBE-vdW functional. This represents an increase of 67 kJ mol−1 as compared to the PW91 functional and agrees with previous theoretical values by Yildirim et al.66 using the same functional (−168 kJ mol−1) but obtained at a lower coverage, and also with the experimental range of adsorption enthalpies (−140–−196 kJ mol−1).80 For hydrogen adsorption, however, the adsorption enthalpy calculated with the optPBE-vdW functional (−102 kJ molH2−1) is 22 kJ mol−1 weaker than with the PW91 functional (−124 kJ molH2−1). Although the value for the optPBE-vdW is closer to the experimentally reported −86 kJ molH2−1,74 the weaker adsorption enthalpy is attributed to an incorrect description of the gas phase H2 molecule with optPBE-vdW; the atomization energies of H2 obtained with the PW91 and optPBE-vdW functionals are, respectively, −1.4 and 18.2 kJ mol−1 larger than the value obtained using the CBS-QB3 method,83,84 which is well-known for its accuracy for gas phase species.85 The product adsorption energies are not compared here since these are adapted in the simulations to make the overall gas phase equilibrium match the thermodynamic gas phase equilibrium as calculated from the NIST database (see section 2.2). In the next section, the hollow dominant path is further evaluated with the optPBE-vdW functional considering surfaces with increasing hydrogen coverage.

3.2. Kinetics as a function of coverage

In this section, first the kinetic and thermodynamic parameters for the reactions along the dominant path are evaluated for different hydrogen coverages. Next, a coverage-dependent microkinetic model is constructed that is used to calculate coverages and reaction rates under actual reaction conditions. The desorption of the intermediates 1,3-cyclohexadiene and cyclohexene is also considered in the model to evaluate their possible competition with the third (reaction c2, Fig. 3) and fifth (reaction f1, Fig. 3) hydrogen addition, respectively.
3.2.1. Kinetics and thermodynamics along the dominant ortho path. The kinetics of the reactions of the dominant ortho path on the hollow site are evaluated at four different hydrogen coverages (θH,initial) coadsorbed with the hydrocarbon intermediate. These four cases, with the coadsorbed hydrogen shown in red, are illustrated in Fig. 6. The first case (a), which is named θH,initial = 0, assumes that no hydrogen is coadsorbed with the hydrocarbon. This relates to a hydrocarbon coverage of θBHi = 0.11 and a coverage θH,initial = 0.11 of fcc hydrogen adsorbed in separate unit cells, which is the same approach as discussed above. The three subsequent cases with increasing hydrogen coverage are (b) the reactant hydrocarbon coadsorbed with a hydrogen coverage θH,initial of 0.11, (c) with θH,initial = 0.44 and (d) with θH,initial = 0.67. These are modeled with 1, 4, and 6 hydrogen atoms coadsorbed in the (3 × 3) unit cell, respectively, as shown in Fig. 6b–d. These hydrogen atoms are adsorbed at the most stable sites for hydrogen coadsorption with benzene as calculated in a previous study on benzene–hydrogen coadsorption.15 These are an fcc hollow site for θH,initial = 0.11, three fcc and a single hcp site for θH,initial = 0.44, and 6 hydrogen atoms adsorbed on fcc sites for θH,initial = 0.67, respectively. For hydrocarbons other than benzene, hydrogen is assumed to occupy the same sites as for benzene–hydrogen coadsorption.
image file: c7cy00962c-f6.tif
Fig. 6 Benzene coadsorbed at a hollow-hcp(0) site together with a hydrogen coverage of (a) θH,initial = 0 (both reactants in different unit cells), (b) θH,initial = 0.11, (c) θH,initial = 0.44 and (d) θH,initial = 0.67 in the same unit cell: side and top views of the (3 × 3) unit cell used (top and middle panel); top view of the surface overlayer at each coverage (bottom panel). Pd atoms in the second layer are colored in grey and adsorbed hydrogen atoms are in red to improve readability.

The increased hydrogen coverage results in weaker adsorption of the intermediates benzene, 1,3-cyclohexadiene, and cyclohexene (see Fig. 7). Only for cyclohexane the adsorption enthalpy is largely independent of H-coverage, which is related to the weak interactions between the physisorbed cyclohexane and the adsorbed hydrogen atoms (for details see Tables S8 and S9, and for adsorption geometries see Fig. S2 of the ESI). The adsorption enthalpy of benzene on a surface pre-covered with θH,initial = 0.67 of hydrogen (−77.5 kJ mol−1) is close to the −58 to −70 kJ mol−1 estimated by Chou and Vannice under reactive conditions on Pd powder,26 obtained from parameter estimation in a kinetic model to experimental benzene hydrogenation results, and is within 5 kJ mol−1 of the estimate of −82 kJ mol−1 of Singh and Vannice on Pd/Al2O3, obtained in the same way.38 To the best of our knowledge, no adsorption enthalpies have been reported on Pd(111) for cyclohexane, cyclohexadiene or cyclohexene. If compared to other metals, the calculated adsorption enthalpy of cyclohexane on Pd (−78 to −84 kJ mol−1) is more exothermic than the experimental values derived from TPD experiments on Pt(111), −59 kJ mol−1,86 and on Ni(111), −55 kJ mol−1.87 For cyclohexene, the range (−66 to −109 kJ mol−1) agrees with the coverage-dependent heats of adsorption on Pt(111) obtained from calorimetric measurements at ∼100 K, ranging from −130 kJ mol−1 at a low coverage to −47 kJ mol−1 at a higher coverage.88 The adsorption entropies of the studied hydrocarbons (see Fig. 7b) show no systematic variation with the hydrogen coverage. There is almost no effect of coverage; only for the tightest bonded molecule at the highest hydrogen coverages (benzene at θH,initial = 0.67) does the coadsorbed hydrogen lead to greater mobility of the adsorbed hydrocarbon, resulting in a less negative adsorption entropy.


image file: c7cy00962c-f7.tif
Fig. 7 DFT adsorption enthalpy (a, kJ mol−1) and entropy (b, J mol−1 K−1) of benzene, hydrogen, 1,3-cyclohexadiene, cyclohexene and cyclohexane as a function of hydrogen coverage on the surface. For hydrogen adsorption, the values correspond to the integral adsorption energy of the adsorption of θH,initial hydrogen on a clean surface, expressed as per H2 molecule (optPBE vdw-DF, direct DFT data).

Also, hydrogen adsorption weakens with increasing hydrogen coverage, but it is less pronounced than for benzene (see Fig. 7) (without co-adsorbed hydrocarbon). The integral hydrogen adsorption enthalpy becomes less exothermic with increasing coverage by about 15 kJ mol−1, in line with the decrease of some 10 kJ mol−1 in the same coverage range as reported previously.73,89,90

The influence of hydrogen coverage on the kinetic and thermodynamic parameters of the six hydrogenation steps of the dominant path has also been evaluated, for the same studied hydrogen coverages θH,initial = 0, 0.11, 0.44 and 0.67. For all the hydrogenation steps, the reaction enthalpy and activation energy decrease with increasing hydrogen coverage (see Fig. 8 for the overview and ESI Tables S8–S11 for the detailed values). This can be explained by the stronger destabilizing effect of the coadsorbed hydrogen atoms on the reactants than on transition states and products: the transition state and product are less tightly bound to the surface than the reactant and therefore interact less with the coadsorbed hydrogen atoms, and furthermore the coverage coadsorbed with the product θH,final is lower than the θH,initial coadsorbed with the reactant. The missing hydrogen will be replenished in the adsorption step after the hydrogen addition step. The effect of hydrogen coverage on the activation and reaction entropies is rather limited and the effect of hydrogen coverage on the pre-exponential factors is below one order of magnitude, which is much smaller than the enthalpic contribution. With the effect on enthalpy and entropy combined, the calculations predict an increase of the equilibrium coefficients with hydrogen coverage for every hydrogenation step (see the right panel in Fig. 8). Equilibrium coefficients increase by 2 to 6 orders of magnitude from the lowest (θH,initial = 0) to the highest hydrogen coverage studied (θH,initial = 0.67).


image file: c7cy00962c-f8.tif
Fig. 8 Reaction enthalpy (left), forward activation energy (middle) and equilibrium coefficients (right) as a function of hydrogen coverage for the six hydrogen addition steps of the dominant path BHi + θH,initial ↔ BHi+1 + θH,final, for the four studied θH,initial hydrogen coverages coadsorbed with the hydrocarbon.
3.2.2. Construction of the coverage-dependent microkinetic model. Based on the obtained data a coverage-dependent microkinetic model can be constructed. The enthalpy profile along the dominant ortho path can be plotted for all the 4 hydrogen coverages θH,initial studied, see Fig. 9. The six hydrogen adsorption steps associated with replenishment of the reacted hydrogen on the surface are explicitly shown between consecutive hydrogenation steps, since these are associated with 6 different adsorption enthalpies. The data presented in the graph pertain to the following DFT calculations: first, benzene and a coverage θH,initial of hydrogen adsorb in one (3 × 3) unit cell, followed by sequential additions of a single hydrogen, and, after each addition step, replenishment of the previously reacted coadsorbed hydrogen by adsorption of a single hydrogen atom. Comparing the graphs at different hydrogen coverages, the hydrogen adsorption clearly becomes weaker with increasing coverage, while at the same time the hydrogenation steps become less endothermic. While the hydrogen adsorption enthalpy amounts to −125 kJ molH2−1 on a clean surface, it decreases to −58 kJ molH2−1 on average for the θH,initial = 0.67 case, corresponding to adsorption in a (3 × 3) unit cell precovered with θH = 0.56 and one hydrocarbon molecule. This range of adsorption enthalpies agrees with the value of −70 kJ molH2−1 estimated by Aben et al.24 on Pd/SiO2 fitting a kinetic model to experimental benzene hydrogenation data, but is stronger than the questionable −5 to −8 kJ mol−1 of Chou and Vannice, which was also determined by model fitting.26
image file: c7cy00962c-f9.tif
Fig. 9 Enthalpy diagram of the dominant path of “hollow” benzene hydrogenation for the 4 hydrogen coverages studied (optPBE-vdW, 450 K): reactants adsorbed in different unit cells (black line), hydrocarbons coadsorbed with θH = 0.11 (dashed black line), 0.44 (grey line) and 0.67 (dashed grey line). Hydrogen adsorption is shown between the hydrogenation steps, also for the first case (θH,initial = 0) for which the hydrogen adsorption enthalpy is actually unique. For this case, there is also an initial hydrogen adsorption step shown since no initial hydrogen is present on the surface.

Furthermore, hydrogen tends to adsorb more strongly with increasing degree of hydrogenation of the hydrocarbon present on the surface by 10 to 30 kJ mol−1 for the hydrogen replenishment step from the first to the last intermediate (see Table S9 for details). This can be explained in terms of the decreased repulsive interactions between adsorbed hydrogen and hydrocarbon: from benzene, which is bound to 3 Pd atoms, over the partially hydrogenated intermediates, which are bound to 1 to 2 Pd atoms, and finally to the fully hydrogenated physisorbed cyclohexane, the hydrocarbon is less strongly attached to the surface, which leads to less steric interaction with hydrogen and, hence, stronger Pt–H bonds. This leads to differences in hydrogen adsorption enthalpy that can amount up to ∼30 kJ mol−1 between the first and the last steps.

The microkinetic model requires a single hydrogen species, in which the thermochemistry and kinetics are independent of the actual neighboring adsorbates. However, currently there are six subsequent hydrogen adsorption steps after each hydrogenation step to replenish the hydrogen surface coverage (θH,final + θBHi + ½H2θH,initial + θBHi). Therefore, the unique hydrogen adsorption enthalpy (and entropy) required in the microkinetic model is obtained by averaging the values for the six consecutive H adsorption steps (see Table S9 of the ESI for details). The thermodynamic consistency of the overall reaction is not affected since the sum of the 6 individual adsorption enthalpies is equal to six times the average adsorption enthalpy. For the θH,initial = 0 case, the calculated hydrogen adsorption energy is already unique and averaging is not required.

To allow a microkinetic simulation in which the thermodynamic and kinetic parameters are a smooth function of coverage, the discrete values at θH,initial = 0, 0.11, 0.44 and 0.67 should be cast in a continuous function. To this end, a linear regression is made to the available data as a function of the total coverage θtotal, i.e. the sum of the hydrocarbon and the hydrogen coverages. The linear function of the kinetic parameters as a function of the total coverage θtotal, as given by eqn (7), is shown in Fig. 10 as an example for the forward activation energies. The total coverage represents the fraction of occupied sites on the surface (θtotal = 1 − θ*), which for the four different coverages studied relates to i) θtotal = 0.11 for the case of θBHi = 0.11 and θH = 0.11 adsorbed in different unit cells (θH,initial = 0), ii) θtotal = 0.22 for θBHi = 0.11 and θH = 0.11 adsorbed in the same unit cell (θH,initial = 0.11), iii) θtotal = 0.55 for θBHi = 0.11 and θH,initial = 0.44, and finally iv) θtotal = 0.78 for θBHi = 0.11 and θH,initial = 0.67. The linear character of the correlation ensures the preservation of thermodynamic consistency; the overall gas phase equilibrium remains independent of the actual coverages on the surface. All the coefficients obtained from the regression are shown in Tables S12 and S13 in the ESI.


image file: c7cy00962c-f10.tif
Fig. 10 Activation energies as a function of coverage for the 6 subsequent hydrogen addition steps in the dominant ortho path (hollow site) and linear correlations.
3.2.3. Simulation of the catalyst performance. The coverage-dependent microkinetic model constructed in the previous section is integrated into a differential reactor model to predict activities under actual reaction conditions. The results of the simulations can be compared to previous experimental observations from the literature.26,28,81,82,91,92

The TOF from this coverage-dependent model is orders of magnitude larger than the TOF predicted based on the low-coverage kinetics. The calculated turnover frequency at 413 K under the experimental conditions used by Chou and Vannice26 is 2.5 × 10−1 s−1, which is in the middle of the available experimental data for a range of Pd catalysts on various supports; Chou and Vannice observed 6 × 10−2–1.1 s−1 for various supported catalysts and 6 × 10−3 s−1 for Pd powder26 (see Table 1 and the triangles in Fig. 11a) (pB,in = 6.7 × 10−2 bar; pH2,in = 0.91 bar, evaluated in a CSTR reactor with W/FB,in = 13.5 kgcat s mol−1 resulting in X < 5%). Also, under other conditions reported by other authors, the model simulations perform well (see Table 1).24,26,28,81,82,91,92 As for the experimental data of Chou and Vannice,26 and also for those of Figueras et al.,81 the simulated TOF is within the reported experimental range. For the data of Moss et al.,91 Fuentes et al.92 and Vannice and Neikam,82 the simulated TOF is higher than the experimentally observed TOFs. For the data of Moss et al. and Fuentes et al., this remains limited to about a factor of 3. For the data of Vannice and Neikam, however, the simulated TOF is a factor of 17 higher than the experimental value. The conditions are, however, very similar to most of the other experiments, and the reported TOF is an outlier in the set of experimental values available. The cited experimental literature is not always clear on the presence of side products such as cyclohexene and 1,3-cyclohexadiene; usually the TOF is derived from the conversion of benzene, without specifying selectivities. The simulated data largely validate this assumption, with all but one of the simulated cyclohexane selectivities above 89%.


image file: c7cy00962c-f11.tif
Fig. 11 (Top) Benzene hydrogenation turnover frequency (TOF) (s−1) as a function of the temperature simulated for the coverage-dependent kinetic model. Experimental TOFs are shown for benzene hydrogenation over the different supported Pd catalysts of Moss et al. (●),91 Chou and Vannice (△),26 and over unsupported Pd powder (▲).26 (Bottom) Simulated fractional surface coverage of all surface intermediates (W/FB,in = 13.5 kgcat s mol−1; pB,in = 6.7 × 10−2 bar; pH2,in = 0.91 bar, which correspond to the conditions of Chou and Vannice and almost those of Moss et al., where pB,in = 2 × 10−2 bar) (kinetic parameters from Tables S12 and S13 of the ESI).
Table 1 Simulated turnover frequencies (TOF in s−1) compared to experimental data from the literature, evaluated at the cited temperatures and inlet pressures. All simulations are performed at a space time of W/FB0 = 13.5 kgcat s mol−1
Catalyst T (K) p B (bar) p H2 (bar) Experimental TOF (s−1) Simulated TOF (s−1) Simulated SCHA
a The number of active sites was reported to be determined from CO chemisorption. b Pd supported on several supports: TiO2, SiO2–Al2O3, Al2O3, SiO2. c No reactant pressures reported in the paper (most likely pH2 ≈ 1 bar and pB ≪ 1 bar). d Reactant pressures were not reported but derived from the description of the experimental set-up and flows. e Pd supported on several supports: Al2O3, SiO2, SiO2–Al2O3, MgO, X and Y faujasite. f Interpolation to 413 K by Figueras et al. from measurements in the range 393–443 K (other supports).
Chou and Vannice26a Pd powder 413 6.7 × 10−2 0.91 6 × 10−3 2.6 × 10−1 89%
Supported Pdb 413 6.7 × 10−2 0.91 6 × 10−2–1.1 2.6 × 10−1 89%
Orozco and Webb28c Pd/SiO2; Pd/Al2O3 480 n/a n/a 4.2 × 10−1; 1.4 × 10−1
Vannice and Neikam82d Pd/Al2O3 423 4.3 × 10−2 1.013 2 × 10−2 3.5 × 10−1 90%
Figueras et al.81 Supported Pde 413f 7 × 10−2 0.9 7 × 10−2–2.8 × 10−1 2.7 × 10−1 89%
Fuentes et al.92 Pd/SiO2; Pd/Al2O3 413 7.5 × 10−2 0.94 2.6 × 10−2; 8.5 × 10−2 2.7 × 10−1 89%
Moss et al.91d Pd/SiO2 373 2 × 10−2 1 1.5 × 10−2–8 × 10−2 1.1 × 10−1 70%


As illustrated in Fig. 11, the temperature at which the simulated turnover frequency reaches a maximum is 475 K. At the same reactant pressures and a slightly larger conversion, Chou and Vannice found the maximum to be located near 495 K (see Table 2).26 Orozco and Webb,28 who performed experiments using a hydrogen flow at atmospheric pressure with intermittent injections of benzene, found maxima at 460 and 480 K, for Pd/SiO2 and Pd/Al2O3, respectively. The coverage-dependent kinetic model presented in this work is clearly able to capture the temperature behavior of the activity, with 475 K in the middle of the range of available experimental data. The only difference is that after the maximum, the simulated TOFs decrease slightly more rapidly than the experimentally observed ones (see Fig. S3 for details).

Table 2 Benzene hydrogen turnover frequencies (TOF in s−1), temperature Tmax at which the maximum TOF is located (K), and partial reaction orders for benzene (m) and hydrogen (n), and comparison to existing experimental data. The simulated results are obtained under the conditions of Chou & Vannice, i.e. T = 413 K; W/FB0 = 13.5 kgcat s mol−1; FB0 = 6.7 × 10−2 mol h−1; FH20 = 0.91 mol h−1, pB0 = 6.7 × 10−2 bar; pH20 = 0.91 bar (kinetic parameters from Tables S12 and S13 of the ESI)
Catalyst TOF (413 K) T max R = kpBmpH2na (T = 413 K)
s−1 K m n
a Partial reaction orders evaluated at constant inlet pressures. b Pd supported on several supports: TiO2, SiO2–Al2O3, Al2O3, SiO2. c No reactant pressures reported by these authors (most likely pH2 ≈ 1 bar and pB ≪ 1 bar).
This work Pd(111) 2.5 × 10−1 475 0.01 1.22
Chou and Vannice26 Pd powder 6 × 10−3 495 0.2 1
Supported Pdb 6 × 10−2–1.1 495 −0.5 to 0.15 ∼1
Orozco and Webb28c Pd/SiO2; Pd/Al2O3 2 × 10−1; 7 × 10−2 460–480


Partial reaction orders in benzene and hydrogen were evaluated under the conditions used by Chou & Vannice and compared to their partial reaction orders obtained on various supported Pd catalyts.26 At 413 K, a zero order in benzene is simulated, which is within the reported range of −0.5 to 0.2 by Chou and Vannice (see Table 2). Also, for the partial reaction order in hydrogen, the calculated value of 1.22 is close to the experimentally observed value of 1 at 413 K. The strong variation of the reaction orders with temperature is reproduced well for hydrogen (see Fig. 12). For benzene, the experimentally observed temperature variation of the reaction order is strongly different for unsupported Pd (Pd powder) and supported Pd, most likely due to the possible involvement of the support in adsorbing benzene. Not surprisingly, the simulated reaction order corresponds best to the data obtained on the unsupported Pd, starting from negative reaction orders at 350 K and rising to about 0.5. Only, the experimentally measured order decreases to about 0.1 above 475 K, the temperature of the maximal TOF, while the calculated value shows a steady increase. This could be related to the simulated TOFs decreasing more rapidly after the maximum TOF than the experimental TOFs, which will increase the simulated partial reaction orders.


image file: c7cy00962c-f12.tif
Fig. 12 Partial reaction orders in benzene (top) and hydrogen (bottom) as a function of temperature, simulation (line) vs. the experimentally measured partial reaction orders of Chou & Vannice26 for Pd powder (△) and TiO2, SiO2 and Al2O3 supported Pd (closed rhombi) (W/FB,in = 3.75 kgcat s mol−1; pB,in = 6.7 × 10−2 bar; pH2,in = 0.80 bar) (kinetic parameters from Tables S12 and S13 of the ESI).

Furthermore, the surface coverages of benzene and hydrogen (see bottom panel of Fig. 11) are in good agreement with models derived based on fitting to experimental data;26,27 benzene is the most abundant surface species (40 to 96%) followed by hydrogen (0.1–8%). All other surface coverages are below 10−3%. As the temperature increases, most coverages decrease, except for a steady increase for hydrogen and the BH intermediate, and a maximum for BH3 and adsorbed cyclohexane (BH6) just before the temperature at which the maximal TOF is observed (475 K).

The available experimental literature on benzene hydrogenation focuses on benzene conversion, but it rarely discusses selectivities, assuming all benzene is converted into cyclohexane. In our simulations, cyclohexane indeed is the main product at all temperatures (see Fig. 13). Above 420 K, less than 10% cyclohexene is formed, though at 350 K the model predicts up to 50%, be it at very low turnover frequencies. The selectivity to 1,3-cyclohexadiene remains below 0.02% for all the temperatures studied. This is in line with the study of Derbentsev et al. on benzene hydrogenation over different noble metals,93 in which no cyclohexadiene was found but 0.1 to 5% of cyclohexene was observed (at 450–460 K). Also for benzene hydrogenation over Pt catalysts, a low but non-negligible selectivity towards cyclohexene is observed experimentally on (111) surfaces (S < 10% at pB = 1.3 kPa, pH2 = 13.3 kPa, T = 380–440 K).35 That the simulated cyclohexene selectivity is somewhat larger than what is experimentally observed can have different reasons, but it is most likely related to the competition between desorption and further hydrogenation of cyclohexene. Cyclohexene desorption is modeled faster than experiment because the initial sticking probability s0 is taken as 1 for all the adsorption steps in this work, to avoid inclusion of data from experimental sources. This is, however, significantly larger than the experimentally observed value for cyclohexene adsorption on the structurally similar Pt(111) (0.21 (ref. 94) to 0.95 (ref. 88)). This can affect the selectivity to cyclohexene, since an initial sticking probability of 1 will overestimate the cyclohexene desorption rate coefficient and hence overpredict the cyclohexene desorption rate with respect to its further hydrogenation towards cyclohexane, leading to a higher cyclohexene selectivity. The actual sticking coefficient can be expected to be less than one, which would lead to a lower cyclohexene desorption rate coefficient and hence a lower cyclohexene selectivity.


image file: c7cy00962c-f13.tif
Fig. 13 Turnover frequency for cyclohexane (full line) and cyclohexene (dashed line) during benzene hydrogenation, with the experimental benzene turnover frequencies of Moss et al.91 and Chou and Vannice26 indicated. The TOF of 1,3-cyclohexadiene is below 10−6 s−1 and therefore not shown (W/FB0 = 13.5 kgcat s mol−1; FB0 = 6.7 × 10−2 mol h−1; FH20 = 0.91 mol h−1, pB0 = 6.7 × 10−2 bar; pH20 = 0.91 bar) (kinetic parameters from Tables S12 and S13 of the ESI).

The reactions with the largest influence on the global rate can also be evaluated using sensitivity analysis for the coverage-dependent model. The second hydrogenation step, from monohydrobenzene (BH) to 1,3-cyclohexadiene (13CHD), has the largest influence on the benzene turnover frequency below 425 K (XSA,c1 ∼ 0.5), while above 425 K the overall rate is most sensitive to the third hydrogenation step (XSA,d1 ∼ 0.5) (see Table S14 in the ESI for details). The same holds for the TOFs for cyclohexane and cyclohexene. The main reasons for the shift in rate determining step with increasing temperature are the coverages of the first and second hydrogenation intermediates: the coverage of intermediate BH increases with temperature, accelerating the second step, while the coverage of BH2 decreases with temperature, decreasing the rate of the third step. This simultaneous change causes the rate-controlling step to evolve from the second to the third step with increasing temperature. Apart from these positive sensitivities, there is a negative sensitivity of the cyclohexene TOF to the rate coefficient of the 5th hydrogenation step, and of the cyclohexadiene TOF to the third hydrogenation step. The obvious reason is that increasing the rate coefficient of these reaction steps favours these intermediates to be further hydrogenated rather than desorbed. Overall, the sensitivity results indicate that there is no single reaction that solely controls the global rate, but that the rate-controlling step changes with the conditions. The sensitivity analysis clearly shows that the influence of coverage is essential for identifying the kinetically most important step(s), which makes identification of the rate determining step based on rate coefficients clearly unreliable. Analysis of the rate coefficients only would indicate the first hydrogenation step to be rate determining, which has a high reaction barrier since the aromaticity of benzene is broken in the first hydrogenation step, but the sensitivity of the global reaction to this step is very low due to the benzene coverage being much higher than the coverage of the other reaction intermediates.

In summary, this work has shown that the construction of an explicit coverage-dependent microkinetic model based on optPBE-vdW DFT calculations at increased hydrogen coverages on the surface allows the prediction of catalytic performances that are comparable to experimental observations in terms of TOF, temperature dependence and partial reaction orders. In this coverage-dependent model, the reactions with the largest sensitivity coefficient depend on the reaction conditions, and these reactions are different from those calculated using the low-coverage kinetics or based on rate coefficients.

4. Conclusions

Periodic density functional theory has been used to evaluate the catalytic hydrogenation of benzene over Pd(111). The 26 reactions paths in the two reaction networks that correspond to the two most stable benzene adsorption sites, namely, bridge and hollow-hcp, can be simplified to a single dominant path starting from benzene adsorbed at a hollow-hcp site, which can explain over 90% of the benzene conversion. The dominant path, which remains dominant regardless of temperature, follows the consecutive hydrogenation of carbon atoms in ortho positions relative to the previously hydrogenated carbon atom.

Evaluating the thermodynamic and kinetic parameters as a function of hydrogen coverage in the unit cell for the reactions along the dominant path shows that the adsorption equilibrium coefficients decrease with hydrogen coverage, while the equilibrium and rate coefficients of the hydrogenation steps increase with coverage. Using these data as a starting point, a coverage-dependent microkinetic model is constructed based on linear correlations of the activation and adsorption energies/entropies as a function of the total coverage. Coupling the microkinetic model in a reactor model allows the evaluation of the kinetics under actual reaction conditions. The TOF predicted with this model is of the same magnitude as experimental observations, in contrast to the microkinetic simulations without coverage dependence. Furthermore, the reaction to which the overall rate is most sensitive does not only change with the conditions but is also different from those obtained with the low-coverage model. The final conclusion of the simulations is that properly accounting for coverage dependence on the kinetics and thermodynamics of catalytic hydrogenation reactions is not only required for an accurate DFT-based prediction of the catalytic activity but also for a correct understanding of the reaction mechanism.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The authors acknowledge financial support from the IFP Energies Nouvelles, the Long Term Structural Methusalem Funding by the Flemish Government, the Fund for Scientific Research Flanders (FWO) and the Interuniversity Attraction Poles Programme – Belgian State – Belgian Science Policy (IUAP). This work was carried out using the STEVIN Supercomputer Infrastructure at Ghent University, funded by Ghent University, the Flemish Supercomputer Center (VSC), the Hercules Foundation and the Flemish Government – department EWI.

Notes and references

  1. J. A. Dumesic, D. F. Rudd, L. M. Aparicio, J. E. Rekoske and A. A. Treviño, The Microkinetics of Heterogeneous Catalysis, American Chemical Society, Washington DC, 1993 Search PubMed .
  2. M. K. Sabbe, M. F. Reyniers and K. Reuter, Catal. Sci. Technol., 2012, 2, 2010–2024 CAS .
  3. Y. K. Park, P. Aghalayam and D. G. Vlachos, J. Phys. Chem. A, 1999, 103, 8101–8107 CrossRef CAS .
  4. O. Lytken, W. Lew and C. T. Campbell, Chem. Soc. Rev., 2008, 37, 2172–2179 RSC .
  5. H. Ihm, H. M. Ajo, J. M. Gottfried, P. Bera and C. T. Campbell, J. Phys. Chem. B, 2004, 108, 14627–14633 CrossRef CAS .
  6. D. H. Mei, M. Neurock and C. M. Smith, J. Catal., 2009, 268, 181–195 CrossRef CAS .
  7. A. C. Lausche, A. J. Medford, T. S. Khan, Y. Xu, T. Bligaard, F. Abild-Pedersen, J. K. Norskov and F. Studt, J. Catal., 2013, 307, 275–282 CrossRef CAS .
  8. M. K. Sabbe, G. Canduela-Rodriguez, M. F. Reyniers and G. B. Marin, J. Catal., 2015, 330, 406–422 CrossRef CAS .
  9. E. Shustorovich and H. Sellers, Surf. Sci. Rep., 1998, 31, 5–119 CrossRef .
  10. A. B. Mhadeshwar and D. G. Vlachos, J. Catal., 2005, 234, 48–63 CrossRef CAS .
  11. M. Maestri and K. Reuter, Angew. Chem., Int. Ed., 2011, 50, 1194–1197 CrossRef CAS PubMed .
  12. S. Singh, S. Li, R. Carrasquillo-Flores, A. C. Alba-Rubio, J. A. Dumesic and M. Mavrikakis, AIChE J., 2014, 60, 1303–1319 CrossRef CAS .
  13. C. Wu, D. J. Schmidt, C. Wolverton and W. F. Schneider, J. Catal., 2012, 286, 88–94 CrossRef CAS .
  14. C. Chizallet, G. Bonnard, E. Krebs, L. Bisson, C. Thomazeau and P. Raybaud, J. Phys. Chem. C, 2011, 115, 12135–12149 CAS .
  15. G. Canduela-Rodriguez, M. K. Sabbe, M. F. Reyniers, J. F. Joly and G. B. Marin, Phys. Chem. Chem. Phys., 2014, 16, 23754–23768 RSC .
  16. A. C. Capleton and L. S. Levy, Chem.-Biol. Interact., 2005, 153, 43–53 CrossRef PubMed .
  17. M. Sanati, B. Harrysson, M. Faghihi, B. Gevert and S. Jaras, in Catalysis, ed. J. J. Spivey, Royal Society of Chemistry, Cambridge, 2002, pp. 1–42 Search PubMed .
  18. H. A. Franco and M. J. Phillips, J. Catal., 1980, 63, 346–354 CrossRef CAS .
  19. M. A. Keane and P. M. Patterson, Ind. Eng. Chem. Res., 1999, 38, 1295–1305 CrossRef CAS .
  20. C. Mirodatos, J. A. Dalmon and G. A. Martin, J. Catal., 1987, 105, 405–415 CrossRef CAS .
  21. R. Z. C. Vanmeerten and J. W. E. Coenen, J. Catal., 1977, 46, 13–24 CrossRef CAS .
  22. R. Z. C. Van Meerten, T. F. M. Degraaf and J. W. E. Coenen, J. Catal., 1977, 46, 1–12 CrossRef CAS .
  23. R. Z. C. Van Meerten, A. C. M. Verhaak and J. W. E. Coenen, J. Catal., 1976, 44, 217–225 CrossRef CAS .
  24. P. C. Aben, J. C. Platteeu and B. Stoutham, Recl. Trav. Chim. Pays-Bas, 1970, 89, 449 CrossRef CAS .
  25. L. P. Lindfors, T. Salmi and S. Smeds, Chem. Eng. Sci., 1993, 48, 3813–3828 CrossRef CAS .
  26. P. Chou and M. A. Vannice, J. Catal., 1987, 107, 129–139 CrossRef CAS .
  27. P. Chou and M. A. Vannice, J. Catal., 1987, 107, 140–153 CrossRef CAS .
  28. J. M. Orozco and G. Webb, Appl. Catal., 1983, 6, 67–84 CrossRef CAS .
  29. A. Palazov, T. Bonev and D. Shopov, React. Kinet. Catal. Lett., 1978, 9, 383–387 CrossRef CAS .
  30. S. D. Lin and M. A. Vannice, J. Catal., 1993, 143, 539–553 CrossRef CAS .
  31. M. Yang, K. C. Chou and G. A. Somorjai, J. Phys. Chem. B, 2003, 107, 5267–5272 CrossRef CAS .
  32. K. M. Bratlie, L. D. Flores and G. A. Somorjai, J. Phys. Chem. B, 2006, 110, 10051–10057 CrossRef CAS PubMed .
  33. K. M. Bratlie, L. D. Flores and G. A. Somorjai, Abstr. Pap. Am. Chem. Soc., 2006, 231, COLL88 Search PubMed .
  34. K. M. Bratlie, C. J. Kliewer and G. A. Somorjai, J. Phys. Chem. B, 2006, 110, 17925–17930 CrossRef CAS PubMed .
  35. K. M. Bratlie, M. O. Montano, L. D. Flores, M. Paajanen and G. A. Somorjai, J. Am. Chem. Soc., 2006, 128, 12810–12816 CrossRef CAS PubMed .
  36. J. P. G. Kehoe and J. B. Butt, J. Appl. Chem. Biotechn., 1972, 22, 23 CrossRef CAS .
  37. M. Xi and B. E. Bent, J. Vac. Sci. Technol., B: Microelectron. Nanometer Struct.--Process., Meas., Phenom., 1992, 10, 2440–2446 CrossRef CAS .
  38. U. K. Singh and M. A. Vannice, AIChE J., 1999, 45, 1059–1071 CrossRef CAS .
  39. J. W. Thybaut, M. Saeys and G. B. Marin, Chem. Eng. J., 2002, 90, 117–129 CrossRef CAS .
  40. S. L. Lu, W. W. Lonergan, J. P. Bosco, S. R. Wang, Y. X. Zhu, Y. C. Xie and J. G. Chen, J. Catal., 2008, 259, 260–268 CrossRef CAS .
  41. M. Saeys, M. F. Reyniers, M. Neurock and G. B. Marin, J. Phys. Chem. B, 2003, 107, 3844–3855 CrossRef CAS .
  42. P. Tetenyi and Z. Paal, Z. Phys. Chem., 1972, 80, 63 CrossRef CAS .
  43. X. C. Su, K. Y. Kung, J. Lahtinen, Y. R. Shen and G. A. Somorjai, J. Mol. Catal. A: Chem., 1999, 141, 9–19 CrossRef CAS .
  44. D. P. Land, C. L. Pettiettehall, R. T. Mciver and J. C. Hemminger, J. Am. Chem. Soc., 1989, 111, 5970–5972 CrossRef CAS .
  45. Y. I. Derbentsev, Z. Paál and P. Tétényi, Z. Phys. Chem., 1972, 80, 51–62 CrossRef CAS .
  46. T. Takahashi, K. Yamashita, T. Kai and I. Fujiyoshi, Can. J. Chem. Eng., 1986, 64, 1008–1013 CrossRef CAS .
  47. M. Saeys, M. F. Reyniers, M. Neurock and G. B. Marin, J. Phys. Chem. B, 2005, 109, 2064–2073 CrossRef CAS PubMed .
  48. F. Mittendorfer and J. Hafner, J. Phys. Chem. B, 2002, 106, 13299–13305 CrossRef CAS .
  49. C. Morin, D. Simon and P. Sautet, Surf. Sci., 2006, 600, 1339–1350 CrossRef CAS .
  50. P. Horiuti and M. Polanyi, Trans. Faraday Soc., 1934, 30, 1164–1172 RSC .
  51. T. Bera, J. W. Thybaut and G. B. Marin, Ind. Eng. Chem. Res., 2011, 50, 12933–12945 CrossRef CAS .
  52. G. Kresse, J. Non-Cryst. Solids, 1995, 193, 222–229 CrossRef .
  53. G. Kresse and J. Furthmuller, Comput. Mater. Sci., 1996, 6, 15–50 CrossRef CAS .
  54. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186 CrossRef CAS .
  55. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269 CrossRef CAS .
  56. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979 CrossRef .
  57. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775 CrossRef CAS .
  58. M. Methfessel and A. T. Paxton, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40, 3616–3621 CrossRef CAS .
  59. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192 CrossRef .
  60. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed .
  61. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6687 CrossRef CAS .
  62. H. Rydberg, B. I. Lundqvist, D. C. Langreth and M. Dion, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 6997–7006 CrossRef CAS .
  63. M. Dion, H. Rydberg, E. Schroder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401 CrossRef CAS PubMed .
  64. J. Klimeš, D. R. Bowler and A. Michaelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 22, 022201 CrossRef PubMed .
  65. J. Klimeš, D. R. Bowler and A. Michaelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 195131 CrossRef .
  66. H. Yildirim, T. Greber and A. Kara, J. Phys. Chem. C, 2013, 117, 20572–20583 CAS .
  67. S. Gautier, S. N. Steinmann, C. Michel, P. Fleurat-Lessard and P. Sautet, Phys. Chem. Chem. Phys., 2015, 17, 28921–28930 RSC .
  68. G. Mills, H. Jonsson and G. K. Schenter, Surf. Sci., 1995, 324, 305–337 CrossRef CAS .
  69. G. Henkelman and H. Jonsson, J. Chem. Phys., 1999, 111, 7010–7022 CrossRef CAS .
  70. National Institute of Standards and Technology Chemistry webbook, Standard Reference Database 69, June 2005 Release, http://webbook.nist.gov. .
  71. G. B. Marin and G. S. Yablonsky, Kinetics of chemical reactions : decoding complexity, Wiley-VCH, Weinheim, 2011 Search PubMed .
  72. J. A. Konvalinka and J. J. F. Scholten, J. Catal., 1977, 48, 374–385 CrossRef CAS .
  73. J. F. Paul and P. Sautet, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 8015–8027 CrossRef CAS .
  74. H. Conrad, G. Ertl and E. E. Latta, Surf. Sci., 1974, 41, 435–446 CrossRef .
  75. C. Morin, D. Simon and P. Sautet, J. Phys. Chem. B, 2004, 108, 5653–5665 CrossRef CAS .
  76. M. Saeys, M. F. Reyniers, G. B. Marin and M. Neurock, J. Phys. Chem. B, 2002, 106, 7489–7498 CrossRef CAS .
  77. F. Mittendorfer, C. Thomazeau, P. Raybaud and H. Toulhoat, J. Phys. Chem. B, 2003, 107, 12287–12295 CrossRef CAS .
  78. M. K. Sabbe, L. Lain, M. F. Reyniers and G. B. Marin, Phys. Chem. Chem. Phys., 2013, 15, 12197–12214 RSC .
  79. W. T. Tysoe, R. M. Ormerod, R. M. Lambert, G. Zgrablich and A. Ramirezcuesta, J. Phys. Chem., 1993, 97, 3365–3370 CrossRef CAS .
  80. A. F. Lee, K. Wilson, R. M. Lambert, A. Goldoni, A. Baraldi and G. Paolucci, J. Phys. Chem. B, 2000, 104, 11729–11733 CrossRef CAS .
  81. F. Figueras, R. Gomez and M. Primet, Adv. Chem. Ser., 1973, 480–489 CrossRef CAS .
  82. M. A. Vannice and W. C. Neikam, J. Catal., 1971, 23, 401–405 CrossRef CAS .
  83. J. A. Montgomery, M. J. Frisch, J. W. Ochterski and G. A. Petersson, J. Chem. Phys., 1999, 110, 2822–2827 CrossRef CAS .
  84. J. A. Montgomery, M. J. Frisch, J. W. Ochterski and G. A. Petersson, J. Chem. Phys., 2000, 112, 6532–6542 CrossRef CAS .
  85. P. D. Paraskevas, M. K. Sabbe, M. F. Reyniers, N. Papayannakos and G. B. Marin, Chem. – Eur. J., 2013, 19, 16431–16452 CrossRef CAS PubMed .
  86. C. Xu, Y. L. Tsai and B. E. Koel, J. Phys. Chem., 1994, 98, 585–593 CrossRef CAS .
  87. P. Zebisch, W. Huber and H. P. Steinruck, Surf. Sci., 1991, 244, 185–196 CrossRef CAS .
  88. O. Lytken, W. Lew, J. J. W. Harris, E. K. Vestergaard, J. M. Gottfried and C. T. Campbell, J. Am. Chem. Soc., 2008, 130, 10247–10257 CrossRef CAS PubMed .
  89. W. Dong, V. Ledentu, P. Sautet, A. Eichler and J. Hafner, Surf. Sci., 1998, 411, 123–136 CrossRef CAS .
  90. O. M. Lovvik and R. A. Olsen, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 10890–10898 CrossRef .
  91. R. L. Moss, D. Pope, B. J. Davis and D. H. Edwards, J. Catal., 1979, 58, 206–219 CrossRef CAS .
  92. S. Fuentes and F. Figueras, J. Chem. Soc., Faraday Trans. 1, 1978, 74, 174–181 RSC .
  93. Y. Derbentsev, P. Tetenyi and Z. Paal, Z. Phys. Chem., 1972, 80, 51–62 CrossRef CAS .
  94. F. C. Henn, A. L. Diaz, M. E. Bussell, M. B. Hugenschmidt, M. E. Domagala and C. T. Campbell, J. Phys. Chem., 1992, 96, 5965–5974 CrossRef CAS .

Footnote

Electronic supplementary information (ESI) available: (i) Calculation of rate and equilibrium coefficients, (ii) procedure for the correction of gas-phase thermodynamics, and thermodynamic & rate parameters and sensitivity coefficients for (iii) the low-coverage kinetic model and (iv) the coverage-dependent kinetic model. See DOI: 10.1039/c7cy00962c

This journal is © The Royal Society of Chemistry 2017