Density Functional Theory Based Micro- and Macro-Kinetic Studies of Ni-Catalyzed Methanol Steam Reforming

: The intrinsic mechanism of Ni-catalyzed methanol steam reforming (MSR) is examined by considering 54 elementary reaction steps involved in MSR over Ni(111). Density functional theory computations and transition state theory analyses are performed on the elementary reaction network. A microkinetic model is constructed by combining the quantum chemical results with a continuous stirring tank reactor model. MSR rates deduced from the microkinetic model agree with the available experimental data. The microkinetic model is used to identify the main reaction pathway, the rate determining step, and the coverages of surface species. An analytical expression of MSR rate is derived based on the dominant reaction pathway and the coverages of surface species. The analytical rate equation is easy to use and should be very helpful for the design and optimization of the operating conditions of MSR.


Introduction
Fuel cell is an attractive environmentally friendly energy conversion technology, but its application is very limited due to the difficulties associated with the transportation and storage of hydrogen [1][2][3]. A promising route for the broad adoption of fuel cells is through the so-called onboard fuel cell technology where H 2 is obtained via real-time steam reforming of hydrocarbon, such as methanol, which is liquid at room temperature and normal pressure [3,4].
Cu is the most commonly used commercial catalyst for the methanol steam reforming (MSR), CH 3 OH + H 2 O CO 2 + 3H 2 . However, Cu-catalyst suffers from some serious drawbacks such as pyrophoricity and catalyst sintering [5][6][7]. To overcome the drawbacks of Cu catalyst, numerous alternative materials have been examined, including noble metal catalysts [8][9][10], Ni-Cu alloy-based catalysts [11][12][13], and Ni-based catalysts [7,[14][15][16]. On one hand, noble metals (Pd, Pt) are known to be the most active catalysts for MSR [8][9][10], but the high prices limit their large-scale industrial applications. On the other hand, Ni-Cu alloys are shown to be better than Cu as the MSR catalyst. For example, Lytkina et al. showed that Ni 0.2 Cu 0.8 /ZrO 2 had the highest activity and best CO 2 selectivity among their Ni-Cu alloy catalysts [11]. Khzouz et al. attributed their observed high MSR activity of Cu-Ni alloy to the higher methanol decomposition activity of Ni than Cu [12]. Together with the observed high CH 3 OH conversion activities of Ni-catalysts [14][15][16][17], and considering their low prices, Ni-based catalysts appear to be promising industrial catalysts for MSR.

Quantum Chemical Results of Elementary Reactions
A total of 15 surface species are considered. The most stable adsorption sites and adsorption energies of the adsorbates are listed in Table 1. As shown in Table 1, H 2 O*, CH 3 OH*, and CO 2 * are only weakly bounded on Ni(111). The diffusion barrier of CO 2 * is found to be about 3.85 kJ/mol, which is less than k B T when T > 200 • C. Therefore, CO 2 * is viewed as moving freely on Ni(111). The rotational barriers of H 2 O* and CH 3 OH* are 1.92 kJ/mol and 0.96 kJ/mol, respectively. Although the adsorption energy of CH 3 O* is high, its rotational barrier is only 3.85 kJ/mol. Therefore, the three surface species of H 2 O*, CH 3 O*, and CH 3 OH* are allowed for 1D-rotational motions. The information about the diffusional and rotational motions of surface species is used in the computation of Q R in Equation (20).  Table 2, where A is a temperature independent pre-exponent constant. As seen in Table 2, the energy barriers of multiple reaction pathways are comparable, preventing a simple determination of the main reaction pathway and the rate-determining step (RDS). Therefore, constructing some micro-kinetic model is required to reveal the main reaction pathway and RDS and make a comparison with experiments.  (1).

Comparison of Microkinetics with Experiment
The rate constant parameters in Table 2 are used in the reaction rate per unit surface area of the species and combined with the ideal CSTR model to simulate the MSR kinetics. First, the microkinetic model is compared with the experiment of Liu et al. [7]. The experimental conditions of steam methanol ratio of 3:1 and a flow rate of 10 cc/min are used in the microkinetic model. Moreover, the total catalyst surface area is assumed to be 0.14 m 2 , an estimation based on the experimental description of 2-3 mg of catalyst sample and catalyst particle size of 12 nm. The theoretical and experimental comparison thus obtained is shown in Figure 1. As seen in Figure 1, the theoretical and experimental methanol conversion rates show similar variation patterns and are in a semi-quantitative agreement. The maximal and average theoretical and experimental differences in the methanol conversion rate are about 15% and 12%, respectively. Notice that the theoretical reaction rates are very sensitive to the computed parameters such as the activation energies; it is quite often to see order(s) of magnitude difference in the theoretical and experimental kinetic results. For example, the reaction rates of Ni-catalyzed methane steam reforming as predicted by two DFT-based microkinetic models are 1 to 3 orders of magnitude below the measured values [23,24]. Even for a kinetic model obtained by fitting experimental data, an order of magnitude difference may be found when applying the model to the experimental data that are not used in the model fitting [25]. Therefore, the theoretical and experimental conversion rates shown in Figure 1 can be viewed as in a very good qualitative agreement.
Catalysts 2020, 10, 349 4 of 11 quantitative agreement. The maximal and average theoretical and experimental differences in the methanol conversion rate are about 15% and 12%, respectively. Notice that the theoretical reaction rates are very sensitive to the computed parameters such as the activation energies; it is quite often to see order(s) of magnitude difference in the theoretical and experimental kinetic results. For example, the reaction rates of Ni-catalyzed methane steam reforming as predicted by two DFT-based microkinetic models are 1 to 3 orders of magnitude below the measured values [23,24]. Even for a kinetic model obtained by fitting experimental data, an order of magnitude difference may be found when applying the model to the experimental data that are not used in the model fitting [25]. Therefore, the theoretical and experimental conversion rates shown in Figure 1 can be viewed as in a very good qualitative agreement.

The Reaction Flux
Because the intermediate energies of multiple pathways are similar with comparable energy barriers, as inferred from Table 2, the flux analysis is employed to discern the main reaction pathway. The results of flux analysis for a representative operating condition is shown in Figure 2. As seen in Figure 2a, nearly all reaction flux flows through the main pathway. The reaction pathway can be understood by considering the data in Table 2 as the following: As the forward activation energy of Reaction (4) is far lower than that of Reaction (2)-(4) contributes to almost all decomposition of methanol after its adsorption on Ni(111). CH3O* is then dehydrogenated in two steps to yield CHO*. CHO* is further dehydrogenated to produce CO* via Reaction (15). The consumption of CHO* via Reaction (16) to form HCOO** is negligible due to its relatively high activation barrier and small prefactor. This is quite different from the case of MSR on Cu/ZnO surface, where the pathway related to the formation of HCOO** is very important [5,26,27]. CO* is then converted into CO2* mostly through Reaction (17), CO*+O*→CO2*, even though the forward activation reactions of Reaction (17) and Reaction (21) are comparable. This is due to the decomposition of COOH* into CO2* encountering a high energy barrier, resulting in a low reaction flux for the pathway of Reaction (21) and Reaction (22), CO*+OH*→COOH*→CO2*+H*.
The forward and reverse reaction fluxes of the main elementary reaction pathway of MSR are shown in Figure 2b. As seen in Figure 2b, while CH3OH+* ⇄ CH3OH* and CO2* ⇄ CO2+* are two equilibrium steps, the remaining reactions are nonequilibrium. The nonequilibrium steps complicate the understanding of MSR kinetics, and a simple, yet reasonably accurate model, is highly desirable.

The Reaction Flux
Because the intermediate energies of multiple pathways are similar with comparable energy barriers, as inferred from Table 2, the flux analysis is employed to discern the main reaction pathway. The results of flux analysis for a representative operating condition is shown in Figure 2. As seen in Figure 2a, nearly all reaction flux flows through the main pathway. The reaction pathway can be understood by considering the data in Table 2 as the following:  Table 2.

Derivation of Analytical Macrokinetic Rate Expression
As seen in Figure 2b, the reverse reaction rate of CH3O*+* ⇄ CH2O*+H* is only 0.7% of the forward reaction. Ignoring the reverse flux only overestimates the net rate of Reaction (13) Table 2.

of 11
As the forward activation energy of Reaction (4) is far lower than that of Reaction (2)-(4) contributes to almost all decomposition of methanol after its adsorption on Ni(111). CH 3 O* is then dehydrogenated in two steps to yield CHO*. CHO* is further dehydrogenated to produce CO* via Reaction (15). The consumption of CHO* via Reaction (16) to form HCOO** is negligible due to its relatively high activation barrier and small prefactor. This is quite different from the case of MSR on Cu/ZnO surface, where the pathway related to the formation of HCOO** is very important [5,26,27]. CO* is then converted into CO 2 * mostly through Reaction (17), CO*+O*→CO 2 *, even though the forward activation reactions of Reaction (17) and Reaction (21) are comparable. This is due to the decomposition of COOH* into CO 2 * encountering a high energy barrier, resulting in a low reaction flux for the pathway of Reaction (21) and Reaction (22), CO*+OH*→COOH*→CO 2 *+H*.
The forward and reverse reaction fluxes of the main elementary reaction pathway of MSR are shown in Figure 2b. As seen in Figure 2b, while CH 3 OH+* CH 3 OH* and CO 2 * CO 2 +* are two equilibrium steps, the remaining reactions are nonequilibrium. The nonequilibrium steps complicate the understanding of MSR kinetics, and a simple, yet reasonably accurate model, is highly desirable.

Derivation of Analytical Macrokinetic Rate Expression
As seen in Figure 2b, the reverse reaction rate of CH 3 O*+* CH 2 O*+H* is only 0.7% of the forward reaction. Ignoring the reverse flux only overestimates the net rate of Reaction (13) by 0.7%, r net 13 = r + 13 − r − 13 = r + 13 − 0.7%r + 13 = 99.3% r + 13 ≈ r + 13 . Moreover, approximating CH 3 OH*+* CH 3 O*+H* as an equilibrium step may only overestimate the CH 3 O* concentration by about (5.1 × 10 −9 -4.7 × 10 −9 )/ 4.7 × 10 −9 = 8.5%. Therefore, it is reasonable to simplify the MSR kinetics as where ⇔ indicates an equilibrium step, while → denotes an irreversible reaction, i.e., zero flux for the reverse reaction. The net reaction rate of methanol is then given by where θ and θ CH 3 O * are respectively the coverages of surface activation sites and CH 3 O*. The equilibrium conditions of Reaction (1) and Reaction (3) yield, where P i (I = CH 3 OH) denotes the partial pressure of the gas species i. Combining Equations (3)- (5) gives The flux analysis shows the H 2 decomposition, Reaction (24), is an equilibrium step and one has Consequently, The coverage of surface activation sites depends on the coverages of surface species. Table 3 shows the coverages of surface species corresponding to the operating parameters of Figure 2. The dominant surface species are CO*, OH*, O*, and H*, with a combined coverage of about 96.8%. The coverage of surface activation sites is 3.2%. All the other surface species cover less than 0.001% of the surface. Although the coverages of surface species vary with the operating condition, the coverages of the minority surface species are found to be very small in all simulated cases. Based on the above observation, one may write, As microkinetics simulations show that Reactions (19), (20), (25), and (26) are always equilibrium steps, we have, Equations (7), (10)-(12) yield, Inserting the coverages of surface species into Equation (9), the coverage of surface activation sites is found to be Combining Equation (16) with Equation (8) yields the reaction rate of MSR, Based on the data in Table 2, the rate parameters in Equation (17) can be calculated as  (17) is obtained by the forward reaction of MSR. The overall MSR rate is found by multiplying Equation (17) with a reaction kinetic factor of F kin = (1− where K MSR is the equilibrium constant of MSR, For normal MSR conditions, F kin ≈ 1 is expected, and the difference between Equations (17) and (18) is negligible.
Equation (17) (Equation (18)) is easy to use. Once the gas compositions and the temperature are known, the reaction rate of MSR can be calculated. Moreover, the coverages of main surface species are also known through Equation (7) and Equations (13)-(16).

DFT Calculations of Elementary Reactions
Quantum chemical computations are carried out with the DFT computational software Vienna Ab-initio Simulation Package (VASP.5.4.1, University of Vienna, Vienna, Austria, 2016) [28][29][30][31]. The projector augmented wave (PAW) method [32,33] is used to describe the electron-ion interaction between core ion and valence electrons. The cutoff energy of wavefunctions of valence electrons is 380 eV. The exchange-correlation energy functional of OptB88-vdW is chosen, as it is the best for describing the van der Waals (vdW) interaction on a metal surface [34]. As the size of an Ni-particle in most applications is in the order of 1 micron, the (111) surface should be the dominant catalytic surface [24]. Therefore, Ni(111) is used to build the surface model that consists of a 3 × 3 unit cell with a three-layer slab and a vacuum region of 10 Å. The Brillouin zone is sampled by a 5 × 5 × 1 k-point Monkhorst-Pack grid.
DFT calculations are used to optimize the structures of MSR intermediates on Ni(111). The adsorption energy and the most stable sites of intermediates is then obtained. The adsorption energy is defined here as: ∆E ads where E x_sur f is the energy of adsorbate and surface, E x is the energy of isolate x molecule/atom, and E sur f is the energy of clean Ni(111) surface. In the geometric optimization, the convergence criteria for electronic computation is 1 × 10 -6 eV/atom, and all geometries were optimized using an energy-based conjugate gradient algorithm until ionic energy was converged below 1 × 10 -5 eV/atom.
In the transition state calculations, considering both efficiency and accuracy, a saddle point is searched in two steps. First, the computationally efficient climbing image nudged elastic band (CL-NEB) method [35,36] is used to locate the minimum energy path and the transition state. In CL-NEB calculations, 6 images are inserted between the initial state and the end state, the electron energy convergence criteria is 1 × 10 -4 eV, and the force is converged to 1 × 10 -7 eV/atom. Second, the CL-NEB transition states are used as the initial structures of the dimer method [37] to find the accurate transition state structures. In the dimer method calculations, the electron and force convergence are increased respectively to 1 × 10 -7 eV/atom and 0.01 eV/Å. Vibrational frequencies of the relevant structures are calculated based on the simple harmonic oscillator approximation. In the frequency calculations, the displacement per degree of freedom is 0.005 Å and the electronic energy is converged to 1 × 10 -8 eV/atom. The vibrational frequencies are used to compute the reaction constants of the elementary reactions as well as to validate the transition states.

Micro-Kinetic Model
The microkinetics of MSR elementary reactions on Ni(111) is constructed based on the TST theory and the results of DFT computations. A continuous-stirred tank reactor (CSTR) is used to simulate the microkinetic reaction process and coupled with a flux analyses [38] to gain understandings of reaction pathways and kinetic behaviors. A CSTR model is used as it is a reasonable reflection of the flow cell used in the Ni-catalyzed MSR experiments [7,39]. Nevertheless, it is noted that the correspondence between an ideal CSTR model and the experimental flow cell is inexact and only semi-quantitative in nature.
In a CSTR model, the mass conservation equation is given by, where f n , W n , and r n are the mass fraction, the molar weight, and the molar production rate per unit surface area of the nth species, respectively, Y n,in is the mass fraction of the nth species in the inlet gas, . m in and . m out are respectively the inlet and outlet mass flow rate, ρ is the mass density of gas in the tank reactor, V is the reactor volume, and A s is the surface area of the catalyst.
The reaction rate per unit surface area of the nth species can be calculated as where C m is the concentration of the mth gas-phase or surface species involved in the reaction, k l,f and k l,r are the forward and reverse reaction rate constants of the lth reaction, respectively, v f l,n is the stoichiometric coefficient of the nth species in the lth reaction, and the superscripts f and r indicate respectively the forward and reverse stoichiometric coefficients.
According to the TST [40], the rate constant in Equation (20) is where k B is the Boltzmann constant, T is the temperature, h is the Plank constant, Q' is the partition function of the transition state excluding the motion along the reaction coordinate, E a is the activation energy, R is the ideal gas constant, and Q R is the partition function of the reactant. Here, Q R is dependent on the motion modes of the reactant. The motion modes of a surface species are related to its translational/rotational energy barrier, E a . For k B T > E a , the vibrational mode of adsorbate in the direction of rotation/translation can be converted into a translational/rotational mode. There are no activation energy barriers for the adsorption of some gas molecules. The rate constant for such non-activated adsorption is calculated as [23] k ads = k B T 2πM (22) where M is the mass of the adsorbed molecule. The reaction rate constant for molecular desorption with no desorption energy barrier is computed by the definition of equilibrium constant k des = k ads K eq (23) where the equilibrium constant, K eq , is the ratio of the partitioning functions of the adsorbed and desorbed molecules [40]. The reaction rate constant in Equation (21) is temperature dependent. The temperature dependent rate constant is fitted to the analytical form Equation (1). The fitting is made for T = 200-500 • C.

Conclusions
A comprehensive set of elementary reaction steps of MSR on Ni(111) are studied by DFT calculations using the OptB88-vdW functional. A microkinetic model is built based on the DFT results and the TST theory. The microkinetic model is in good agreement with the available experimental data. Microkinetic simulations of MSR on Ni(111) show clearly: (1) CH 3 OH → CH 3 OH * → CH 3 O * → CH 2 O * → CHO * → CO * → CO 2 * → CO 2 is the main reaction pathway, (2) CH 3 O*→CH 2 O*+H* is the RDS, and (3) O*, CO*, OH*, and H* are the only surface species with non-negligible surface coverages.
Based on the information revealed by the microkinetic simulations, an analytical macrokinetic rate equation is derived. The analytical rate equation of MSR is not only easy to use, but also very helpful for guiding the design and optimization of MSR operating conditions.