First-Principles Microkinetic Modeling of Methane Oxidation over Pd

Research Article pubs.acs.org/acscatalysis

First-Principles Microkinetic Modeling of Methane Oxidation over Pd(100) and Pd(111) Mikkel Jørgensen* and Henrik Grönbeck* Department of Physics and Competence Centre for Catalysis, Chalmers University of Technology, 412 58 Göteborg, Sweden

Downloaded via UNIV OF CALIFORNIA SANTA BARBARA on January 6, 2019 at 18:21:20 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The intrinsic activity of Pd(100) and Pd(111) for methane oxidation is investigated by Density Functional Theory (DFT)-based mean-field microkinetic modeling. The model includes 32 reaction steps, and the calculated turnover frequencies together with reaction orders compare favorably with experimental data. On both surfaces, the reaction proceeds via complete dehydrogenation of methane to elemental carbon followed by different mechanisms for carbon oxidization. Pd(100) is found to be more active than Pd(111) at temperatures from 400 to 1000 K. For both surfaces, the reaction order in methane approaches unity with increasing temperature. The reaction order in water is positive at low temperatures owing to water-promoted carbon oxidation. Methane dissociation is the main rate-controlling step for Pd(111), whereas formation of COH and CO is also controlling the rate over Pd(100). The present work uncovers the detailed reaction mechanisms for complete methane oxidation over palladium, which can be used in catalyst design to target the rate-controlling steps. KEYWORDS: heterogeneous catalysis, microkinetic modeling, DFT, methane oxidation, palladium, Pd(100), Pd(111)

■

INTRODUCTION Methane is known to be a potent greenhouse gas with a global warming potential exceeding that of CO2.1,2 It is therefore important to limit methane emissions which is generally done by complete catalytic oxidation to CO2 and H2O. Pd and Pt are known to be efficient catalysts for Complete Methane Oxidation (CMO), where Pd is superior in oxygen-rich environments and Pt has the highest activity under oxygenpoor conditions.3 As palladium is oxidized during oxidizing conditions, it has been difficult to establish the most active phase of the catalyst during methane oxidation. Some studies have reported metallic Pd as the most active phase,4−6 whereas other results indicate that the high activity should be attributed to PdO.7,8 It has also been argued that a mixed Pd/PdO phase is advantageous for a high catalytic activity. The reason for the debate is that metallic Pd is unstable under relevant O2 pressures and temperatures.2,7,9 The experimental accessibility to metallic Pd surfaces or nanoparticles is furthermore obscured by the tendency of Pd to form surface oxides.9,10 In fact, it has been reported that the saturation coverage of oxygen is 0.25 monolayers (ML) and that surface oxides are thermodynamically preferred at higher coverages.9 A (√5 × √5) R 27° oxide-overlayer has been determined for Pd(100),9 whereas a (√6 × √6) surface oxide is grown on Pd(111).10 Dissociative methane adsorption to methyl and hydrogen is generally assumed to be the rate-determining step in CMO over metallic Pd.1,11−15 This is a reasonable assumption at high temperatures owing to low surface coverages and facile oxygen adsorption. The assumption could probably be questioned at © 2016 American Chemical Society

lower temperatures where the presence of different adsorbates could alter the mechanism. Previous studies have usually postulated a reaction mechanism where methane is dissociated via a pyrolytic pathway to hydrogen and elemental carbon.6,16−18 The reaction proceeds thereafter along one branch to CO2 by sequential addition of oxygen to carbon, and along another branch by oxidation of hydrogen to H2O. An alternative mechanism has been suggested, which is based on kinetic isotope experiments and Density Functional Theory (DFT) calculations,15,18 involving methane dissociation over chemisorbed oxygen to form OH+CH3 or alternatively OH +CH3O. It should be noted that the assumed oxygen coverage in refs 15,18 was up to 0.985 ML which is substantially higher than the experimentally observed saturation coverage.9,10 A recent study of methane oxidation over Pt(111) indicate the importance of intermediates such as CH2OH, COH, CHO, CHOH, and OCOH.19 This is an observation that could also have a bearing on the reaction mechanism on metallic Pd. CMO over PdO has been studied experimentally4,5,13,20−24 as well as theoretically,15,25−27 and it has been suggested that PdO(101) is the most active surface of the oxide. A recent theoretical study on methane oxidation on PdO(101) shows that the reaction mechanism is complex and sensitively dependent on reaction conditions.25 For example, it was found that CO2 is formed via CH2OH from which hydrogen is Received: June 21, 2016 Revised: August 25, 2016 Published: August 31, 2016 6730

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis

species at some conditions.37 The implications of this finding for the reaction kinetics, however, require further investigation. Calculations of reaction barriers are done by use of the climbing image Nudged Elastic Band (NEB)38 method from the VTST tools39 using at least five intermediate images and forces converged to a maximum of 0.05 eV/Å. Initial interpolations between images are done using the Image Dependent Pair Potential method.40 Vibrational analysis of the considered transition states (TS) confirm that the structures are true saddle points having one imaginary frequency parallel to the reaction coordinate. The reaction barriers are evaluated as the energy difference between the TS and the adsorbates in separate cells. Microkinetic Modeling. The reaction kinetics are treated by mean field microkinetic modeling. This involves solving a system of coupled ordinary differential equations describing the fractional surface coverages. Each species that can be present on the surface contributes with one equation to the system, which has the generic form.41

abstracted and further carbon oxidation proceeds via a Mars− van Krevelen mechanism. Given the recent detailed investigations on CMO for PdO(101), it becomes important to advance the understanding of methane oxidation also over metallic palladium. In particular, there is a need to explore alternative reaction paths both in the decomposition of CH4 and in the formation of CO2 and H2O. Despite the rapid development of experimental in situ techniques, it is still a considerable challenge to control the phase of the catalyst and determine the reaction mechanisms. Instead, this is conveniently done by computational screening of different plausible reaction mechanisms. In the present work, we combine first-principles calculations and kinetic modeling to investigate CMO on Pd(100) and Pd(111) including 32 reaction steps. The model allows us to study complete oxidation without substantially constraining the possible reaction mechanism. To the best of our knowledge, this is the first study of this type for complete methane oxidation over palladium. A key result of the present work is that it reveals competing reaction mechanisms and that the dominant pathway depends on the reaction conditions.

dθi = dt

■

COMPUTATIONAL METHODS DFT Method. Total energies and vibrational frequencies are calculated using DFT28,29 by use of the Plane-Wave Vienna Abinito simulation package (VASP).30−33 The interaction between the valence electrons and the cores are described by the Projector-Augmented Wave (PAW) scheme.34 The number of electrons treated variationally are Pd(10), O(6), C(4) and H(1). The calculations are performed within the Generalized Gradient Approximation (GGA) using the PBE35 exchange− correlation functional. Oxygen chemisorption at a 0.25 ML coverage on Pd(111) is used for tests of convergence with respect to plane-wave cutoff, k-point sampling, and vacuum width. The oxygen adsorption energy is found to be converged within 0.05 eV with a planewave cutoff of 450 eV, a 6 × 6 × 1 k-points grid, and a 12 Å vacuum. The surface slabs are represented by four atomic layers which is sufficient for converged surface energies. The surfaces are constructed from the theoretical lattice constant which is calculated to be 3.956 Å using a (12 × 12 × 12) k-point grid. The setup of atoms, structural optimizations, and vibrational calculations are performed using the Atomistic Simulation Environment (ASE).36 Structures are considered to be optimized when all forces are below 0.01 eV/Å. In these calculations, the two bottom Pd layers are kept fixed to emulate a truncated bulk surface. The energies of gas phase molecules are calculated in a (30 Å × 30 Å × 30 Å) cell with the proper spin states (triplet for O2 and singlet for CO, CO2, H2O, and CH4). Adsorbed species are considered to be singlets. The vibrational modes are calculated from the optimized geometries using the harmonic approximation and finite differences (two-point difference with steps of 0.01 Å) in ASE.36 For adsorbates, all Pd atoms are fixed during the vibrational analysis. The vibrational energies are used to correct all adsorbate and transition-state energies for zero-point motion. Due to the numerical uncertainty in the calculation of low frequency modes, all wavenumbers lower than 100 cm−1 are set to 100 cm−1. The vibrational frequencies are furthermore used to evaluate the entropy of adsorbed species and the vibrational entropy contribution to gas-phase molecules. It has been suggested recently that the harmonic approximation fails to estimate the total entropy of surface

∑ rk(θ)nki k

(1)

Where θi is the fractional coverage of species i, rk is the rate of reaction k that depends on the coverages θ. nki is the number of molecules of species i consumed in reaction k. The fractional coverage is defined by considering each Pd atom as a site. Thus, θi is the ratio of species i to the number of Pd atoms. The system of differential equations is numerically integrated until steady state is reached using SciPy42 with the VODE solver for stiff problems. Rate constants are calculated using harmonic transition-state theory43 by the following expression: kTST =

kBT Z ‡ h Z

(2)

Where kB is Boltzmann’s constant, T is the temperature, h is Planck’s constant, Z‡ is the partition function of the transition state, and Z is the partition function of the initial state. The partition functions in the transition states and of the adsorbates are approximated by treating each degree of freedom as a frustrated vibration. For gas-phase molecules, we estimate the partition functions by statistical thermodynamics, see the Supporting Information. For the oxygen desorption energy, we include a coveragedependent linear O−O repulsion using values of ref 17. This is important as oxygen turns out to be an abundant surface species. A repulsive dependence is also included for the C−C interaction, which is elaborated on in the Supporting Information. Furthermore, we have included OH−OH interactions following the method of ref 17, although this does not affect the model due to low OH coverages. The methane dissociation barrier is also given an oxygen coverage dependence. This coverage dependence is calculated with coverages from 1/9 to 1/4, where the low coverage values are calculated in a (3 × 3) unit cell, see Supporting Information. We choose to include coverage dependences for oxygen and carbon due to their high coverages on the two surfaces. The repulsive adsorbate−adsorbate interactions are affecting the turnover frequencies (TOF) mainly in the low-temperature region where the coverages are significant. A sensitivity analysis of the TOFs with respect to the strengths of the lateral interactions are presented in the Supporting Information. 6731

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis

Table 1. Considered Reaction Steps with the Forward (Ef) and Backward (Eb) Energy Barriers (eV) without Zero-Point Correctionsa no.

equation

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30) (31) (32)

CH4(g) + 2*↔ CH3*+H* CH4(g) + O*+*↔ CH3*+OH* CH*3 + *↔ CH*2 +H* CH3*+O*↔ CH2OH*+* CH2* + *↔ CH*+H* CH*2 +O*↔ CH2O*+* CH*2 +OH*↔ CH2OH*+* CH*+*↔ C*+H* O2(g) + *↔ O2* C*+O*↔ CO*+* C*+OH*↔COH*+* CO*+O*↔ CO2(g) + 2* CO*+OH*↔ OCOH*+* CO*+H*↔COH*+* CO*→CO(g) + * O*+H*↔ OH*+* OH*+H*↔ H2O*+* OH*+OH*↔ H2O*+O* H2O*↔ H2O(g) + * CH2OH*+*↔CH2O*+H* CH2OH*+*↔ CHOH*+H* CH2O*+*↔ CHO*+H* CHO* + * ↔ CO * + H* CHO* + O * ↔ CO* + OH* CHO* + * ↔ CH* + O* COH* + O* ↔OCOH* + * OCOH* ↔ CO2(g) + H* CHOH* + * ↔ COH* + H* CHOH* + 2* ↔CO* + 2H* 2H* → H2(g) + 2* O2* + *↔ 2O* 2O* → O2(g) + 2*

E(100) f

A(100) f

0.79 1.10 0.66 1.03 0.17 1.11 1.07 0.65 0.0 1.81 1.39 0.86 0.72 1.62 1.86 0.77 0.76 0.32 0.28 0.66 0.90 0.29 0.33 0.0 2.43 0.0 0.77 0.77 0.46 0.99 0.13 2.57

× × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×

2.1 5.0 1.4 4.8 4.3 2.1 3.8 7.0 1.9 4.2 4.2 1.5 3.8 9.9 1.0 2.1 2.8 3.4 2.5 7.3 1.1 2.0 4.7 8.7 2.6 1.4 2.8 1.2 1.1 9.8 9.5 2.6

2

10 101 1013 1012 1013 1012 1013 1013 10−1 1012 1013 1013 1012 1011 1013 1011 1012 1012 1013 1013 1014 1013 1013 1012 1013 1012 1014 1014 1014 1012 1012 1013

E(100) b

A(100) b

0.52 1.26 0.57 1.59 0.83 1.76 1.28 1.00 1.48 3.02 1.42 1.01 0.20 0.87

5.1 × 10 2.9 × 1013 4.6 × 1011 8.9 × 1013 3.6 × 1011 2.1 × 1012 8.9 × 1013 7.4 × 1011 1.5 × 1013 4.8 × 1012 9.0 × 1012 7.9 × 10° 1.8 × 1013 4.5 × 1013 11

× × × × × × × × × × × × × ×

1013 1013 1011 101 1011 1011 1011 1011 1012 1012 1013 10−1 1011 109

1.20 0.95 0.08 0.0 0.67 1.23 0.99 1.19 1.40 1.74 0.76 1.9 1.27 1.70

5.1 2.6 1.3 7.5 1.3 6.9 7.8 1.2 5.5 5.5 3.7 1.2 1.4 3.0

1.10

9.7 × 1011

E(111) f

A(111) f

0.99 1.41 0.98 1.65 0.32 1.46 1.00 1.26 0.0 1.70 0.95 1.60 1.09 1.96 1.99 1.16 0.92 0.0 0.23 0.88 0.87 0.39 0.41 0.0 2.06 0.0 0.90 0.53 0.80 1.52 0.55 2.80

× × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×

4.8 3.3 2.7 4.1 3.7 8.3 2.9 3.6 5.9 5.6 1.4 3.2 8.4 2.6 1.0 1.0 6.7 2.4 1.5 1.8 3.1 1.5 8.5 4.1 8.4 1.4 7.8 1.1 2.2 6.2 1.3 9.1

2

10 102 1013 1013 1013 1012 1012 1014 10−1 1012 1013 1013 1012 1013 1013 1014 1013 1013 1013 1014 1014 1014 1013 1013 1012 1012 1014 1013 1014 1015 1013 1013

E(111) b

A(111) b

0.86 1.21 0.95 2.04 0.86 1.94 1.50 0.96 0.77 3.97 2.24 1.52 0.59 0.90

4.5 × 1013 3.3 × 1013 1.8 × 1013 2.0 × 1014 3.9 × 1013 2.0 × 1013 2.1 × 1013 9.2 × 1013 1.5 × 1013 1.7 × 1012 9.1 × 1012 2.7 × 10° 3.4 × 1013 6.3 × 1013 × × × × × × × × × × × × × ×

1014 1013 1012 101 1013 1013 1013 1013 1012 1013 1012 101 1012 1013

1.09 1.22 0.47 0.0 0.95 1.23 1.14 1.70 1.27 1.37 0.59 1.39 1.22 2.54

1.1 2.7 8.9 4.4 5.6 6.6 2.6 1.6 8.0 2.1 2.6 1.5 1.2 1.0

2.44

7.4 × 1013

a The pre-exponential factors (1/s) for forward and backward reactions are denoted Af and Ab, respectively. The pre-exponential factors are temperature dependent and here given for a temperature of 700 K. The pre-exponential factors for adsorption are given in [Pa−1 s−1].

energy arises from lateral interactions in the considered surface cells. For the complete reaction, we calculate an energy change of −7.58 eV, including zero-point corrections. This value could be compared with the experimental gas-phase data of −8.36 eV.44 We find an entropy change at standard conditions of −0.02 meV/K, and the corresponding experimental value is −0.07 meV/K.44

Reaction barriers for nonactivated processes are set to 0.1 eV to avoid numerical issues. The entropic change of barrierless reactions are calculated by taking the average of the entropy in the final and initial state. The estimate of the entropy change is crude; however, the stiffness of the kinetic equations and especially the slow rate for methane adsorption makes the results insensitive to this choice. The kinetic model is ensured to be thermodynamically consistent by comparing the Gibbs free energy change for the overall reaction with that of the surface reactions. This implies that the following equation should be fulfilled:

∑ νi ln(K i) = − i

ΔG kBT

■

RESULTS Potential Energy Landscapes. The considered reaction network with calculated activation energies and pre-exponential factors are reported in Table 1. The investigated reactions are based on suggestions in the literature,17,18,25,45,46 as well as an iterative process where the kinetic modeling provided ideas of alternative reaction steps. Specifically, we included all steps that were considered in ref 17. These are the stepwise methane dehydrogenation in steps 1, 3, 5 and 8, and the oxidation of carbon to CO2 in steps 10 and 12. The water formation via two distinct pathways (steps 16−19) are also taken from ref 17. We included methane dissociation over an oxygen covered surface to form OH in step 2, as previously suggested in ref 47. The steps involving CH2OH, CH2O, and CHO were inspired by the work of ref 25, where methane oxidation was modeled over PdO(101). To our knowledge, the other steps have not been

(3)

Where ΔG is the Gibbs free energy change of the overall reaction, T is the temperature, Ki is the equilibrium constant of reaction i, and νi is the stoichiometric number of reaction i. νi is defined as the number of times reaction i has to proceed to complete one full catalytic cycle. The thermodynamic consistency is evaluated along different reaction paths, and a maximal entropic deviation between the gas-phase and surface reaction occurs at 1000 K and has the value 0.6 eV, which is small at this high temperature. In energy, we find a maximal deviation of 0.20 eV, which is also small. The slight deviation in 6732

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis

where hydrogen is directed toward the surface. We also included an alternative reaction pathway for CO2 formation via CH2OH → CHOH → COH continuing to CO + H or OCOH. This path has been reported to be of importance in methanol dehydrogenation over Pd(100) and Pd(111).45,46 We obtained similar barriers as in refs 45,46,49. Water formation was considered via two mechanisms, namely, step 17, where H and OH react, and step 18, where two OH react to form water. Step 18 is found to be energetically preferred, which is consistent with an earlier DFT study.17 The calculated pre-exponential factors of surface intermediates are generally in the range of 1011 to 1015 s−1, whereas for adsorption, they are in the range 10−1 to 102 Pa−1s−1. The preexponential factors for methane dissociation are consistent with earlier theoretical studies.14,15,17 Experimentally, the preexponential factor has been measured to be 102 Pa−1 s−1 for palladium particles supported on zirconia. A similar preexponential factor has been reported for methane dissociation over Pt.47 Oxygen adsorption was modeled by a molecular precursor (step 9 and step 31), whereas oxygen desorption was modeled in one step (step 32). The calculated pre-exponential factor for oxygen adsorption (step 9) is lower than for methane dissociation. If instead, collision theory and an experimental sticking-coefficient50 of 0.45 is used to estimate the preexponential for O2 adsorption, we obtain ∼102 Pa−1 s−1. This estimate is higher than the value in Table 1 and demonstrates the difficulty in estimating entropy losses in adsorption processes. Fortunately, the reaction rate is not sensitive to this parameter thanks to a rapid oxygen adsorption rate compared to the other relevant reactions. CO-desorption (step 15) has a pre-exponential factor of about 1013 s−1, which is in good agreement with thermal desorption spectroscopy measurements51 that report 1013.5 s−1. For this reaction step, we decided to only include the forward reaction as readsorption is unlikely at low CO pressures. Microkinetic Modeling. Having the energy landscape for methane oxidation over Pd(100) and Pd(111), we turn to the microkinetic modeling. The simulations are performed at oxygen excess conditions similar to the experimental studies of methane oxidation on Pd surfaces.20,21 The methane and oxygen pressures are, thus, set to 0.61 and 3.06 mbar respectively, and the H2O and CO2 pressures are set to zero. The results are not very sensitive to the absolute pressures and the methane/oxygen ratio as long as the oxygen excess condition is preserved. Turnover Frequency and Surface Coverages. The calculated turnover frequency is shown in Figure 2a. The Pd(100) surface is predicted to have the highest activity, where the difference in the activity of the two surfaces is clearest at high temperatures. The TOFs above 900 K are within 1 order of magnitude in agreement with experimental data for Pd(100) and Pd(111).20,21 The TOFs on both surfaces are 1−2 orders of magnitude lower than a previous theoretical study of complete methane oxidation on metallic Pd(100) and Pd(111).17 The difference can mainly be attributed to the assumption in ref 17 that methane dissociative adsorption is the rate-determining step. The microkinetic model gives access to the surface coverages, and the most abundant surface species are reported in Figure 2b,c. The surfaces are mainly covered by O, C, and CO, and the coverage of oxygen ranges from 0.7 to 0.20 ML on both

considered previously in the context of complete methane oxidation. The calculated adsorption energies and reaction barriers are in agreement with previous studies where such exist.17 All structures are included in the Supporting Information. The energy landscape along the CO2 branch for the main reaction pathway at high temperatures is given in Figure 1. The

Figure 1. Reaction energy landscape on Pd(100), top, and Pd(111), bottom, following the carbon atom along the high-temperature pathway. The energies are given with respect to CH4(g) + O2(g) without zero-point corrections. The reaction barriers are indicated. Atom color code: H (white), C (black), O (red), and Pd (blue).

reaction is endothermic until CH forms on the surface. An interesting feature of both surfaces is that the highest barrier is calculated for the oxidation of elemental carbon (step 10) and that the dehydrogenation of methane and hydrocarbyl groups proceeds with fairly low barriers. Even if the energy barrier is low for the methane dissociation (step 1), the reaction could be slow due the large entropy loss upon dissociative adsorption. Elemental carbon is particularly strongly bonded to Pd(100), which has consequences for the reaction kinetics. The COoxidation reaction (step 12) is also identified as a potential kinetic bottleneck at low temperatures due to the rather high barriers of 0.86 and 1.60 eV on Pd(100)and Pd(111), respectively. These values are consistent with a previous DFT study,48 where the barriers were reported to be 0.80 eV on Pd(100) and 1.53 eV on Pd(111). The high barrier for elemental carbon oxidation (step 10) is not compensated by a large pre-exponential factor. This fact led us to introduce an alternative reaction, namely, C + OH → COH (step 11), which turned out to have a lower barrier and higher pre-exponential factor. Further oxidation of COH leads to OCOH which might be kinetically important owing to the high oxygen abundance. The hydrogen abstraction from OCOH to form CO2 has a barrier of 0.77 and 0.90 eV for Pd(100) and Pd(111), respectively. The barrier originates mainly from OCOH isomerization into the TS configuration 6733

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis

methane oxidation over the two surfaces are reported in Figure 3. The apparent activation energy varies considerably with temperature for Pd(100) ranging from 0.84 to 1.20 eV, whereas the value for Pd(111) varies between to 0.88 and 1.13 eV.

Figure 3. Apparent activation energies as a function of temperature of methane oxidation on Pd(100) and Pd(111). The simulations are performed with a methane pressure of 0.61 mbar, an oxygen pressure of 3.6 mbar, and zero water and carbon-dioxide pressure. Figure 2. (a) Natural logarithm of the turnover frequency of methane oxidation over Pd(100) and Pd(111) as a function of temperature. The coverage of O, C, CO, and empty sites on Pd(100) (b) and Pd(111) (c). The simulations are performed with a methane pressure of 0.61 mbar, an oxygen pressure of 3.6 mbar and zero water and carbon-dioxide pressure.

In order to further analyze the difference between the two surfaces at high temperatures, we consider a simple analytical model for the reaction. In particular, only the reaction path discussed in connection to Figure 1 is analyzed with methane dissociation as an irreversible rate-determining step, and all other steps assumed to be in quasi-equilibrium. By doing this, we derive an expression for the apparent activation energy (Eact,app):

surfaces. On Pd(100), carbon is present at the entire temperature range and amounts to almost one-third coverage at 450 K. The reason for the high coverage is the high barriers for the elemental carbon oxidation (step 10) and COH formation from C and OH (step 11). It should be noted that the approach to steady state is slow at low temperatures, in particular, for Pd(100). The steady state over Pd(100) is reached after 106 s starting with an initially clean surface. An integration time of 2 h yields a carbon coverage of 0.002 at 400 K, which should be compared with the steady-state value of 0.27. The slow buildup of carbon may make high carbon coverages difficult to observe experimentally. However, the potentially high coverage together with reasonably low barriers for subsurface diffusion52 suggest that Pd-carbide formation is plausible with carbonaceous adsorbates, which indeed has been observed experimentally.53 On Pd(111), CO is present below 500 K, which agrees with temperature-programmed desorption experiments.54 The carbon coverage is lower than on Pd(100), which is reasonable given that the barrier for COH formation is relatively low on Pd(111). On Pd(111), carbon is present below 700 K, and the steady state is reached after 108 s. An integration time of 2 h at 400 K yields a carbon coverage of zero in this case, which again suggests that it might be difficult to observe adsorbed carbon experimentally. Apparent Activation Energies. A common way to analyze experimentally a reaction is to determine the apparent activation energy. This is a convoluted measure that is determined by the dominant reaction path and may change with reaction conditions. The apparent activation energies for

Eact,app = E1f −

pO K (ΔH9 + ΔH31) 2

1+

pO K 2

(4)

Here, E1f is the methane dissociation barrier, pO2 is the oxygen pressure, K is the equilibrium constant between the forward reaction of step 9 and step 32, which are in equilibrium due to the high temperature. ΔH9 is the heat of adsorption for molecular oxygen adsorption and ΔH31 is the enthalpy change of step 31. The inclusion of the molecular precursor state for O2 has a negligible effect and if this state is omitted, K can be substituted with the equilibrium constant for dissociative adsorption and ΔH9+ΔH31 with the adsorption energy of 2O. The expression for the apparent activation energy shows that Eact,app reflects the oxygen chemisorption energies and barriers for methane dissociation at high temperatures. Using the zero-point corrected values at 1000 K (K = 5 × 10−5 Pa−1, ΔH9 = −1.46 eV, ΔH31 = −1.04 eV and E1f = 0.63 eV) for the Pd(100) surface, yields Eact,app = 0.90 eV, which is close to the value of the numerical simulation which is 0.85 eV. With the corresponding values for Pd(111) (K = 2 × 10−5 Pa−1, ΔH9 = −0.76 eV, ΔH31 = −1.92 eV and E1f = 0.82), the Eact,app is calculated to be 1.01 eV. This should be compared to the value for the numerical simulation which is 0.97 eV. The calculated apparent activation energies are slightly lower than published experimental values for a Pd-foil20−22 that ranged from 1.30 to 1.45 eV. We speculate that the discrepancy might be due to the facile formation of the (√5 × √5) R 27° 6734

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis

The reaction order in water is calculated to be 0.5 at low temperatures on both surfaces. This is attributed to water assisted carbon oxidation. At low temperatures, CO2 is formed via COH, and the abundance of this species is promoted in the presence of water trough the C + OH reaction (step 11). In experiments, it has previously been shown that the presence of water enhances CO-oxidation over Pt(111) at low temperatures.55 At high temperatures the reaction order in water is calculated to be zero which is consistent with the results in refs 20,21. The positive reaction order in water is a characteristic feature of metallic palladium, as oxidized palladium is known to have a clear negative reaction order22,25 due to OH-poisoning. The reaction order in oxygen is calculated to be slightly negative over the entire temperature range for both surfaces ranging from −0.01 to −0.1. The negative sign is a consequence of oxygen blocking sites that are of importance to the dominant reaction path at the given conditions. The reaction order in the CO2 pressure was found to be 0 over the entire temperature range. Preferred Reaction Mechanisms. Having discussed the general behavior of the reaction, we turn to the details in the reaction mechanisms. The preferred reaction mechanisms have been analyzed by inspection of the net rates for the different elementary reactions at steady state. In the limit of high temperature (above 800 K), the reaction follows the scheme outlined in Figure 1. Methane in this path is sequentially dehydrogenated to elemental carbon and hydrogen. CO2 is formed by the C + O and CO + O reactions, whereas water is formed from O + H and OH + OH. The situation is more complicated at lower temperatures where we find that the reaction mainly proceeds via two alternative mechanisms, see Figure 5. The mechanism begins with methane dehydrogenation to C, in similarity to the high temperature mechanism. In mechanism 1, COH is formed from C + OH. In the next step, COH is oxidized by O to OCOH. OCOH is thereafter decomposed to CO + OH which in the

surface oxide over which the barrier for methane oxidation has been calculated to be 1.34 eV.13 Reaction Orders. Experimentally, the reaction orders in methane, oxygen, and water have been determined at high temperatures to probe the intrinsic properties of the metal phase.20,21 The reaction order in methane was measured20 to be 0.9 ± 0.2 and 0.7 ± 0.1, for Pd(100) and Pd(111) respectively. For oxygen, the reaction order was measured to be slightly negative, whereas it was determined to be close to zero for water.20,21 We have analyzed the reaction orders, and the results for methane and water are shown as a function of temperature in Figure 4. On Pd(100), the reaction order in methane pressure approaches unity slowly as the temperature is increased, whereas the approach to unity is more rapid over Pd(111).

Figure 4. Reaction orders in methane and water pressures as a function of temperature over Pd(100) and Pd(111). The simulations are performed with a methane pressure of 0.61 mbar, an oxygen pressure of 3.6 mbar, and zero water and carbon-dioxide pressure.

Figure 5. Two alternative preferred reaction mechanisms at temperatures below 800 K. The simulations are performed with a methane pressure of 0.61 mbar, an oxygen pressure of 3.6 mbar, and zero water and carbon-dioxide pressure. Atom color code: H (white), C (black), O (red), and Pd (blue). 6735

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis

Figure 6. Degree of rate-control analysis for Pd(100) and Pd(111). The simulations are performed with a methane pressure of 0.61 mbar, an oxygen pressure of 3.6 mbar, and zero water and carbon-dioxide pressure.

temperatures, other steps are limiting the rate, such as elemental carbon oxidation, COH formation, and OCOH dehydrogenation. We have additionally identified CO desorption to be inhibiting on Pd(111) at high temperatures due to the comparatively slow CO-oxidation. We find that the most abundant surface species are O, C, and CO. It is well-established that Pd will be oxygen covered in oxidizing environments and that CO will be present on the surface until oxidation light-off. The presence of elemental carbon during the reaction is rationalized by the high barrier for oxidation of elemental carbon, and the possibility to form Pdcarbide has been discussed previously.52,53 The effect of water on the reaction has been frequently discussed in the literature.20−22,25 On metallic Pd, a reaction order between −0.1 and 0.1 has been reported,20,21 whereas other studies22 show a negative reaction order of −0.9. In ref 22, the highly negative order was suggested to be due to PdO formation. A negative reaction order in water over PdO is consistent with a previous DFT-based kinetic model for methane oxidation over PdO(101).25 Here we have shown that the negative order in water pressure appears only on PdO or oxide-overlayers and that the effect of water is negligible on metallic Pd at experimentally relevant temperatures. A positive reaction order in oxygen is also connected to methane oxidation over a PdO phase as the reaction in this case proceeds via a Mars−van Krevelen mechanism consuming lattice oxygen, and adding extra loosely chemisorped oxygen promotes the reaction.25 Over the metal surface, we find that the reaction order in oxygen is slightly negative due to site blocking. The kinetic behavior of the reaction is readily explained by the degree of rate-control analysis. At low temperatures, the reaction has a positive reaction order in water, which in the degree of rate-control analysis is signaled by the high χ for C + OH conversion to COH, where OH originates from dissociated water. Similarly, the reaction order in methane follows closely the temperature dependence of the χ for dissociative methane adsorption. The buildup of an appreciable carbon coverage on Pd(100) is visible in the degree of rate-control analysis trough significant χ values for the oxidation reactions of elemental

last step is oxidized to CO2(g). In mechanism 2, OCOH is directly dehydrogenated to CO2(g). In the intermediate temperature regime, the different pathways are competing on both surfaces. However, the competition is more pronounced on Pd(100) and extends over a temperature interval of 300 K, from 500 to 800 K. This interval is considerably narrower (50 K) on Pd(111), where the high temperature mechanism is dominant already at 550 K. Steps Controlling the Rate. The degree of rate control56 (χi) is a measure of to what extent an elementary reaction step (i) controls the overall rate. We perform this analysis by increasing the forward and backward rate constants of each reaction step by a small amount, and we calculate the linear response in the TOF. This is done while simultaneously keeping the rate of all other reaction steps fixed. The result of this analysis is shown in Figure 6. At low temperatures, Pd(100) has three rate-controlling steps, namely, COH formation through C + OH, CO formation through C + O, and methane dissociation. The reason for having multiple steps controlling the rate is that both the lowand the high-temperature mechanisms are active. However, the low-temperature mechanism has a higher degree of rate control below 500 K. From 500 to 600 K, the CO formation becomes increasingly important. At higher temperatures, methane dissociation (step 1) becomes the sole rate-controlling step, which is in agreement with the experimental literature.1,11,12,15 On Pd(111) below 500 K, the rate-controlling steps are COH formation (step 11) and methane dissociation (step 1). OCOH formation (step 13) is inhibiting the reaction by consuming OH and thereby slowing the formation of COH. The OCOH dehydrogenation to CO2 (step 27) is too slow to compensate for the loss of activity via the dominant pathway. Above 500 K the situation is simpler, where methane dissociation is the rate-determining step. The analysis additionally shows that with increasing temperature, CO desorption is a slightly inhibiting step as it desorbs prior to oxidation. Discussion. In multiple studies, methane dissociation has been assumed to be the rate-determining step in complete methane oxidation.14,15,17 The present work shows that this assumption is valid only at high temperatures. At low 6736

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis

companies AB Volvo, ECAPS AB, Haldor Topsøe A/S, Scania CV AB, Volvo Car Corporation AB, and Wärtsilä Finland Oy.

carbons at low temperatures. The apparent activation energy and the TOF are also connected to the degree of rate control as Eapp,act, to a first approximation, is the sum of the barriers for the reaction steps weighted by the corresponding χ value.

■

■

(1) Gélin, P.; Primet, M. Appl. Catal., B 2002, 39, 1−37. (2) Hoflund, G. B.; Li, Z.; Epling, W. S.; Göbel, T.; Schneider, P.; Hahn, H. React. Kinet. Catal. Lett. 2000, 70, 97−103. (3) Burch, R.; Loader, P.; Urbano, F. Catal. Today 1996, 27, 243− 248. (4) Hicks, R. F.; Qi, H.; Young, M. L.; Lee, R. G. J. Catal. 1990, 122, 280−294. (5) Hicks, R. F.; Qi, H.; Young, M. L.; Lee, R. G. J. Catal. 1990, 122, 295−306. (6) Oh, S. H.; Mitchell, P. J.; Siewert, R. M. J. Catal. 1991, 132, 287− 301. (7) Farrauto, R. J.; Hobson, M. C.; Kennelly, T.; Waterman, E. M. Appl. Catal., A 1992, 81, 227−237. (8) Baldwin, T. R.; Burch, R. Appl. Catal. 1990, 66, 359−381. (9) Lundgren, E.; Gustafson, J.; Mikkelsen, A.; Andersen, J. N.; Stierle, A.; Dosch, H.; Todorova, M.; Rogal, J.; Reuter, K.; Scheffler, M. Phys. Rev. Lett. 2004, 92, 046101. (10) Ketteler, G.; Ogletree, D. F.; Bluhm, H.; Liu, H.; Hebenstreit, E. L. D.; Salmeron, M. J. Am. Chem. Soc. 2005, 127, 18269−18273. (11) Ciuparu, D.; Lyubovsky, M. R.; Altman, E.; Pfefferle, L. D.; Datye, A. Catal. Rev.: Sci. Eng. 2002, 44, 593−649. (12) Li, Z.; Hoflund, G. B. J. Nat. Gas Chem. 2003, 12, 153−160. (13) Martin, N. M.; Van den Bossche, M.; Hellman, A.; Grönbeck, H.; Hakanoglu, C.; Gustafson, J.; Blomberg, S.; Johansson, N.; Liu, Z.; Axnanda, S.; Weaver, J. F.; Lundgren, E. ACS Catal. 2014, 4, 3330− 3334. (14) Hellman, A.; Resta, A.; Martin, N. M.; Gustafson, J.; Trinchero, A.; Carlsson, P.-A.; Balmes, O.; Felici, R.; van Rijn, R.; Frenken, J. W. M.; Andersen, J. N.; Lundgren, E.; Grönbeck, H. J. Phys. Chem. Lett. 2012, 3, 678−682. (15) Chin, Y.; Buda, C.; Neurock, M.; Iglesia, E. J. Am. Chem. Soc. 2013, 135, 15425−15442. (16) Hickman, D. A.; Schmidt, L. D. Science 1993, 259, 343−346. (17) Trinchero, A.; Hellman, A.; Grönbeck, H. Surf. Sci. 2013, 616, 206−213. (18) Chin, Y.; Iglesia, E. J. Phys. Chem. C 2011, 115, 17845−17855. (19) Qi, W.; Ran, J.; Wang, R.; Du, X.; Shi, J.; Ran, M. Comput. Mater. Sci. 2016, 111, 430−442. (20) Zhu, G.; Han, J.; Zemlyanov, D. Y.; Ribeiro, F. H. J. Am. Chem. Soc. 2004, 126, 9896−9897. (21) Zhu, G.; Han, J.; Zemlyanov, D. Y.; Ribeiro, F. H. J. Phys. Chem. B 2005, 109, 2331−2337. (22) Monteiro, R.; Zemlyanov, D.; Storey, J. M.; Ribeiro, F. H. J. Catal. 2001, 199, 291−301. (23) van Giezen, J. C.; van den Berg, F. R.; Kleinen, J. L.; van Dillen, A. J.; Geus, J. W. Catal. Today 1999, 47, 287−293. (24) Fujimoto, K.; Ribeiro, F. H.; Avalos-Borja, M.; Iglesia, E. J. Catal. 1998, 179, 431−442. (25) Van den Bossche, M.; Grönbeck, H. J. Am. Chem. Soc. 2015, 137, 12035−12044. (26) Dianat, A.; Seriani, N.; Ciacchi, L. C.; Bobeth, M.; Cuniberti, G. Chem. Phys. 2014, 443, 53−60. (27) Lv, C.-Q.; Ling, K.-C.; Wang, G.-C. J. Chem. Phys. 2009, 131, 144704. (28) Hohenberg, P.; Kohn, W. Phys. Rev. 1964, 136, B864−B871. (29) Kohn, W.; Sham, L. J. Phys. Rev. 1965, 140, A1133−A1138. (30) Kresse, G.; Furthmuller, J. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (31) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15−50. (32) Kresse, G.; Hafner, J. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561. (33) Kresse, G.; Joubert, D. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758−1775. (34) Blöchl, P. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979.

CONCLUSIONS We have presented DFT-based microkinetic modeling of complete methane oxidation over Pd(100) and Pd(111). The kinetic simulations included 32 reaction steps resulting in numerous competing reaction paths. In the limit of high temperatures, we find that the reaction follows a mechanism where methane is sequentially dehydrogenated to elemental carbon and hydrogen. CO2 is thereafter formed by the C + O and CO + O reactions, whereas water is formed from O + H and OH + OH. At low temperatures, we find that the situation is more complicated, and the oxidation of elemental carbon proceeds via mechanisms including COH and OCOH. Inclusion of these mechanisms turns out to be crucial as these steps are found to control the rate at low temperature. The reaction order in methane is close to one for Pd(111) at temperatures above 600 K, whereas the reaction order approaches unity significantly slower on Pd(100). The reaction order in water is found to be positive at low temperatures, which is owing to water-promoted carbon oxidation. The calculated turnover frequencies and reaction orders are in good agreement with experiments performed over metallic palladium. The most active phase of palladium during methane oxidation has been a long-standing debate in the literature. The difficulty to disentangle the rates for Pd and PdO arises as the oxide phase is thermodynamically preferred in oxygen excess. The present work determines the intrinsic kinetic behavior of metallic palladium, which is important in order to understand and identify the role of different Pd phases. The results provide information that can be utilized in future design of catalysts by the possibility to target rate-controlling steps.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acscatal.6b01752. Details on calculations of partition functions and coverage dependence of activation energies; sensitivity of the TOF to the lateral repulsions (PDF) Adsorbate and transition-state structures (images and .xyz files) together with vibrational wavenumbers (ZIP)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support is acknowledged from the Swedish Research Council and Chalmers Areas of Advance Nanoscience and Nanotechnology and Transport. The calculations were performed at C3SE (Göteborg) and PDC (Stockholm) via a SNIC grant. The Competence Centre for Catalysis (KCK) is hosted by Chalmers University of Technology and is financially supported by the Swedish Energy Agency and the member 6737

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738

Research Article

ACS Catalysis (35) Perdew, J.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (36) Bahn, S. R.; Jacobsen, K. W. Comput. Sci. Eng. 2002, 4, 56−66. (37) Sprowl, L. H.; Campbell, C. T.; Á rnadóttir, L. J. Phys. Chem. C 2016, 120, 9719−9731. (38) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. J. Chem. Phys. 2000, 113, 9901−9904. (39) Transition State Tools for VASP. http://theory.cm.utexas.edu/ vtsttools/ (downloaded/accessed November 2015). (40) Smidstrup, S.; Pedersen, A.; Stokbro, H.; Jónsson, K. J. Chem. Phys. 2014, 140, 214106. (41) Chorkendorff, I.; Niemantsverdriet, J. W. Concepts of Modern Catalysis and Kinetics, 2nd ed.; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, 2007; pp 52−57. (42) Jones, E.; Oliphant, T.; Peterson, P. et al. SciPy, open source scientific tools for Python, 2001; http://www.scipy.org/ (online; accessed February 24, 2016). (43) Eyring, H. Chem. Rev. 1935, 17, 65−77. (44) NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P. J., Mallard, W. G., Eds.; http://webbook.nist. gov; National Institute of Standards and Technology: Gaithersburg, MD, 2016. (45) Jiang, R.; Guo, W.; Li, M.; Lu, X.; Yuan, J.; Shan, H. Phys. Chem. Chem. Phys. 2010, 12, 7794−7803. (46) Jiang, R.; Guo, W.; Li, M.; Fu, D.; Shan, H. J. Phys. Chem. C 2009, 113, 4188−4197. (47) Chin, Y.; Buda, C.; Neurock, M.; Iglesia, E. J. Am. Chem. Soc. 2011, 133, 15958−15978. (48) Vogel, D.; Spiel, D.; Suchorski, Y.; Trinchero, A.; Schlögl, R.; Grönbeck, H.; Rupprechter, G. Angew. Chem., Int. Ed. 2012, 51, 10041−10044. See Supporting Information. (49) Lin, S.; Ma, J.; Ye, X.; Xie, D.; Guo, H. J. Phys. Chem. C 2013, 117, 14667−14676. (50) Peter, M.; Adamovsky, S.; Flores Camacho, J. M.; Schauermann, S. Faraday Discuss. 2013, 162, 341−354. (51) Guo, X.; Yates, J. T. J. Chem. Phys. 1989, 90, 6761−6766. (52) Seriani, N.; Mittendorfer, F.; Kresse, G. J. Chem. Phys. 2010, 132, 024711. (53) Ziemecki, S. B.; Jones, G. A.; Swartzfager, D. G.; Harlow, R. L.; Faber, J., Jr J. Am. Chem. Soc. 1985, 107, 4547−4548. (54) Szanyi, J.; Kuhn, W. K.; Goodman, W. J. Vac. Sci. Technol., A 1993, 11, 1969−1974. (55) Bergeld, J.; Kasemo, B.; Chakarov, D. V. Surf. Sci. 2001, 495, L815−L820. (56) Campbell, C. T. Top. Catal. 1994, 1, 353−366.

6738

DOI: 10.1021/acscatal.6b01752 ACS Catal. 2016, 6, 6730−6738