CO₂ activation dominating the dry reforming of methane catalyzed by Rh(111) based on multiscale modelling

Abstract

The dry reforming of methane (DRM) converts two greenhouse gases (CH₄ and CO₂) to syngas (CO/H₂). Rh-based catalysts are among the most active DRM catalysts, but they still need to be fully understood at the atomic level. In this work, we evaluated the Rh(111)-catalyzed DRM via periodic density functional theory and kinetic Monte Carlo (kMC) simulations, accounting for lateral interactions. The kinetic model consisted of 38 elementary reactions, including adsorption, desorption, and surface chemical reactions. The reaction network considered both the formation of the DRM products and the competitive reverse water-gas shift reaction. kMC simulations indicated direct CO₂ activation takes place, yielding CO* and O*. The CH oxidation path (CH* + O*) was the preferred route to obtain the second CO molecule, and the water formation minimally affected the final H₂/CO ratio. The catalytic system displayed Arrhenius behavior at different temperatures with an apparent activation energy of 53 kJ mol⁻¹. The degree of rate control analysis identified CO₂ activation as the dominant step in Rh(111)-catalyzed DRM, with no evidence of catalyst deactivation. Our study underscores the utility of multiscale modeling for a comprehensive understanding of heterogeneous catalysts from a bottom-up approach.

---

1. Introduction

There is a strong correlation between global warming rise with atmospheric CH₄ and CO₂ levels due to anthropogenic atmospheric emissions.^1,2 From ca. 1950, fossil fuels burning became the dominant source of these anthropogenic emissions, continuously increasing, given the growing energy demand. Moreover, even with renewable energies containing carbon, avoiding greenhouse gas production is not possible.³ Within this context, converting two greenhouse effect gases, CO₂ and CH₄, into valuable syngas (H₂/CO) makes dry reforming of methane (DRM) a highly attractive process both from an academic perspective and for industrial applications to reduce the emission of greenhouse gases.⁴

The main advantage of producing syngas by this route is the formation of lower H₂/CO ratios, typically influenced by the simultaneous occurrence of the reverse water-gas shift (RWGS) reaction, which is preferred for using it in Fischer–Tropsch synthesis to value-added chemicals.^5–7 Nevertheless, DRM faces the challenge caused by the deactivation due to the coke formation on the catalyst's surface.^8–10 Ni-based catalysts are the most convenient for industrial applications due to their low cost and high catalytic activity.^11,12 However, the high carbon depositions on Ni during the reaction^13,14 have suggested group VIII metals as catalysts since they are catalytically active and selective in the DRM.¹⁵ Specifically, Rh-based catalysts exhibited shallow carbon formation. Other studies involving the development of bimetallic Ni–Rh catalysts suggested an improved catalytic activity compared to a pure Ni-based catalyst.¹⁶

Overall, improved DRM catalysts need to be more efficient and coke-resistant, and heterogeneous catalysis plays a prominent role in developing materials capable of controlling the reaction's rate, selectivity, and conversion. In parallel, computational modeling is an alternative to experiments to gather knowledge highly challenging to obtain solely from an empirical approach. However, developing adequate computational models is not straightforward due to heterogeneous catalysis's high complexity and multi-scale nature.^17–19 The elementary processes occur on active catalytic surfaces. However, the structure and composition of active sites might change due to local variations in temperature and concentrations determined by several factors in the reactor scale. Multiscale modeling aims to tackle this challenge.^17,20

The ab initio multiscale modeling approach starts by building up atomic scale-catalytic models and exploring the energetics of elementary steps based on electronic structure calculations. Then, the results are inputs to microkinetic modeling (MKM) approaches. MKM simulations can be performed using a mean-field approximation (MFA)²¹ or kinetic Monte Carlo (kMC) modeling. MFA has a low computational cost but might lead to inaccurate predictions.²⁰ In contrast, kMC simulations can describe the system dynamics accurately since they can directly account for the effect of lateral interactions on surface coverage and reaction rates.^22–26 However, kMC modeling is computationally more demanding than the MFA. In 2011, M. Stamatakis and D. Vlachos²⁴ developed Zacros,²³ a software to run kMC simulations based on the graph-theoretical model that considers the multidentate nature of adsorbates.²⁴ It also incorporates the cluster expansion Hamiltonians to account for the long-range lateral interactions accurately. In other words, the multi-scale approaches allow a more in-depth interpretation of the targeted reactions.

To assess the suitability of noble metals as candidate heterogeneous catalysts for the DRM reaction, we present a multiscale analysis on the Rh(111) surface by combining periodic DFT calculations with kMC simulations following our previous work with the Ru(0001) facet.²⁷ Although the use of noble metals as catalysts is limited by their elevated cost and low availability, the present study provides hints about how they work at the atomic level and which are the critical elementary steps, besides the origin of their high resistance to coke formation, paired with their significant stability and catalytic activity.

First, we use DFT calculations to obtain the energetics of the DRM's main reaction paths catalyzed by the Rh(111) surface. We then perform kMC simulations coupled with a cluster expansion Hamiltonian implemented in Zacros to include lateral interactions between adsorbates. Arising from these simulations and further analysis, we found the main factors governing the kinetics of the overall reaction, which cannot be obtained by solely relying on DFT calculations. While the thermodynamics of the reaction is well understood, knowledge about the reaction mechanism and kinetics is still highly controversial. Therefore, we also used the concept of degree of rate control (DRC) initially proposed by Campbell to identify the rate-determining step (RDS) explained later in this work. Identifying the RDS and a sensitivity analysis can provide valuable ideas for designing novel catalysts.

2. Computational model and method

The computational methods in this work can be divided into density functional theory calculations and kinetic Monte Carlo simulations.

2.1. DFT calculations

All calculations in the present work have been performed at the DFT level using a plane-wave basis set with a pseudopotential projector augmented wave (PAW) method as implemented in the Vienna ab initio simulation package (VASP) code.^28–30 All geometry optimizations were conducted with the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA)³¹ to describe the exchange and correlation effects. The BEEF-vdW functional,³² including van der Waals corrections, was used in our previous work to investigate the overall mechanism of DRM catalyzed by Ru(0001).²⁷ The results showed the BEEF-vdW functional leads to the same energetic trends as those observed with PBE, but at a higher computational cost. Therefore, the energetics and geometries were evaluated only by employing the PBE functional in the current work. Finally, the plane-wave-basis set was truncated at a kinetic energy of 400 eV, and the Brillouin zone was sampled with a 4 × 4 × 1 Monkhorst–Pack k-point mesh.

The Rh(111) surface (ESI† Fig. S1) was modeled by a periodic 3 × 3 four-layer slab with a vacuum separation of 15 Å in the direction perpendicular to the surface to avoid the adsorbates' interactions with their periodic images. All geometry optimization calculations were performed considering a low-coverage of 1/9 monolayer (ML) adsorbates. The bottom layer of the slab was kept fixed at its bulk positions during the geometry optimizations. In contrast, all Rh atoms were fixed for the frequency calculations, and only the adsorbate atoms could move. The energy of isolated molecules was determined at the Γ-point by placing each species in a cubic box with dimensions 15 Å × 15 Å × 15 Å. Each elementary reaction's transition states (TS) were located using the climbing image nudged elastic band method (CI-NEB)³³ and characterized by appropriate frequency analyses, exhibiting only one imaginary frequency. All geometry configurations were allowed to relax using an energy-based conjugate-gradient algorithm until all the forces were smaller than 0.01 eV Å⁻¹. In contrast, a quasi-Newton algorithm was used for the transition-state refinement.

Gibbs energy profiles (see ESI†) were constructed considering no interaction between adsorbates, i.e., considering just one adsorbate per unit cell. Therefore, the Gibbs energies of co-adsorbed species were calculated as the sum of the individual adsorption-free energies. All the Gibbs energies shown in the profiles are relative to the sum of the Gibbs energies of the initial reactants: CH₄ and CO₂.

The main elementary steps of the DRM mechanism were evaluated based on the previously postulated mechanisms in the literature.^34,35 Thermodynamic corrections were calculated using statistical thermodynamics³⁶ at 973.15 K to build the Gibbs energy profile. The Gibbs energy of the gas phase species was calculated by considering translational, rotational, and vibrational contributions to the partition function. Weakly adsorbed species, such as CH₄, CO₂, and H₂, were modeled as bi-dimensional gases, assuming they keep the gas-phase species' full rotational and vibrational modes. Only the vibrational contributions were considered for the remaining chemisorbed species, including the transition state structures. Thermal effects on the metallic surface were not considered.

2.2. kMC simulations

In this work, we employed the graph-theoretical kMC methodology coupled with cluster expansion Hamiltonians and Brønsted–Evans–Polanyi (BEP) relations for the adlayer energetics^23,24,37 as implemented in the software package Zacros (version 2.0).^24,38 This approach has been used before in our previous work,²⁷ thus, here we just briefly summarize their main features.

A kMC simulation shows the temporal evolution of the number and location of adsorbates bound onto a catalyst surface, which is represented as a lattice whose nodes mimic the active sites on which the elementary events occur.²⁰ These elementary events are treated as independent Poisson processes. Thereby, the evolution of the catalytic system is described by the Markovian master equation (MME):³⁹

Considering that the state-to-state transitions (x′ → x) are Markovian,⁴⁰i.e., there is no connection between the state visited before the current state and any state the system will visit next, allows to assign to each elementary event a rate constant k_x′→x. This rate constant represents the probability, per unit of time, that the transition from a state x to a different state x′ occurs. This probability is calculated via TST:where h denotes Planck's constant, k_B is Boltzmann's constant, T is the temperature, and Q^≠ and Q_x (dimensionless) are the quasi-partition functions corresponding to the transition and initial states, respectively.²⁰E^≠_x→x′ is the zero-point corrected energy barrier calculated at the zero-coverage limit from DFT calculations. Hence, k should be multiplied by a statistical factor l that accounts for the several equivalent ways to achieve such a transition state-due to the presence of identical atoms.⁴¹ Finally, during a kMC simulation, the system passes through many states. The collection of all possible such states defines the state space, Ω, to which x and x′ belong, and thus, a kMC trajectory is simply a random walk on Ω.^19,42

In the current work, we defined two kinetically different types of elementary reactions: adsorption/desorption and surface reactions. For a surface Langmuir–Hinshelwood type reaction, the constant rate can be calculated with eqn (2). As the adsorbed species in a surface reaction have only vibrational components, only the vibration quasi-partition function (Q_vib) must be calculated, whereas for gas-phase species, translational, rotational, and vibrational degrees of freedom must be included. The quasi-partition functions can be calculated from standard statistical mechanical expressions,³⁶ using the information from quantum chemistry calculations, i.e., at the DFT level. All energies discussed in this main text include zero-point-energy (ZPE) corrections, and the reference point is the first vibrational energy level (ν = 0). Therefore, the quasi-partition function is expressed without the ZPE contribution (ESI† note 1).

The rates for adsorption of a gas with mass m_X at a given temperature, T, and partial pressure p_X can be calculated through the Hertz–Knudsen equation as:

where A_st corresponds to the effective area of the catalytic site(s) and m_X is the mass of the adsorbing molecule. A_st was calculated as the total area divided by the number of active sites on the surface. The rate of desorption processes can be determined from TST assuming an early 2D gas-like transition state:All adsorption/desorption processes for products of the reaction were considered non-activated and were calculated using the set of eqn (3) and (4), while the methane and carbon dioxide adsorption were considered activated adsorption and were calculated via:The reverse process of the activated adsorption is calculated analogously to eqn (2). Within the kMC simulations, the Rh(111) surface is represented by a two-dimensional 5 × 5 hexagonal periodic lattice-featuring top, hcp, fcc, and bridge adsorption sites – as shown in Fig. 1. In our previous work for the same reaction catalyzed by Ru(0001),²⁷ we found that this unit cell size was enough, as the results did not change upon increasing the lattice size further. Therefore, that additional validation is not addressed here. For all considered reaction intermediates (e.g., OH, CH₃, H…), only one adsorption site was considered to simplify the kMC simulations since they already involve 36 elementary reactions and four site types. Pairwise lateral–lateral interactions considered here were C–C*, CH–CH*, O–O*, C–O*, and C–CH*, since the coverage of their intermediates is significant during the kMC simulations. The coverage of other intermediates remained low during all kMC simulations. Thus, the interaction energies between other adsorbates were not considered.

Fig. 1 This lattice was used to represent the Rh(111) surface in the kMC simulations, depicting the four different types of sites (top, bridge (brg), fcc, and hcp). The dashed line denotes the simulation box, and the thick black line the unit cell.

The range of operating temperatures were chosen according to the experimental conditions employed.¹⁶ All reaction rates (k) are expressed in s⁻¹.

A kinetic Monte Carlo model is constructed for the DRM reaction on the Rh(111) surface based on 19 reversible elementary processes for the molecular mechanism of the DRM, which are described in ESI† Table S1. They include adsorption and desorption of gas-phase species and surface reactions. Several elementary steps considered here have rate constants significantly higher than other elementary steps, resulting in processes that occur at quite different timescales. The fastest processes quickly reach quasi-equilibrium. An unmodified kMC simulation would spend most of the computational time executing uninteresting fast processes rather than the slow relevant reaction events, which are the most relevant ones in estimating catalytic activity. Therefore, their rate constants are scaled to save simulation time and sample the slow relevant events.³⁸ In the present study, we enabled stiffness scaling for all elementary steps to address the separation time scale between the different possible processes, as implemented in Zacros, similarly to our previous work²⁷ (ESI† note 2).

3. Results and discussion

3.1. Energetics

We considered two possible routes for forming the DRM's main and secondary products from the competitive RWGS reaction. All of them share the CO₂ dissociation as an initial step. The CO₂ can dissociate directly into CO* + O* with an energy barrier of 47 kJ mol⁻¹ or can dissociate, assisted by H*, in two steps. The CO₂ activation assisted by H* proceeds first through the COOH* formation and then through the C–O bond cleavage leading to CO* and OH*. The CO₂ reaction with H* produces the COOH* intermediate via an energy barrier of 84 kJ mol⁻¹. The COOH* formation is considered to follow an Eley–Rideal mechanism^43,44 since the CO₂ molecule is practically physisorbed on Rh(111) with an adsorption energy of just −2 kJ mol⁻¹. Besides, CO₂ is weakly chemisorbed (−6 kJ mol⁻¹), when adsorbed on the surface via a bent configuration. This result agrees with an experimental study of the CO₂ dissociation on Rh(111), which suggested that the direct dissociative adsorption for CO₂ is easily carried out at temperatures above 300 K and 1 atm.⁴⁵ Then, the high-energy COOH* intermediate is dissociated into CO* + OH* with a previous slight rotation of the hydroxyl group from pointing down to the surface to pointing away from the surface with an energy barrier (119 kJ mol⁻¹) higher than the direct CO₂ dissociation by 72 kJ mol⁻¹. Thus, direct CO₂ dissociation is the most favorable path to cleave CO_2(g). The first C–H bond activation occurs via an energy barrier of 68 kJ mol⁻¹. After that, the methane cracking is easily obtained through energy barriers lower than 29 kJ mol⁻¹ until CH* forms. The last C–H bond activation from methane (CH* → C* + H*) has the highest energy barrier along the cracking process, being 24 kJ mol⁻¹ higher than the first C–H activation. The C* species can be oxidized by the O* released in the direct CO₂ cleavage to produce CO* through 167 kJ mol⁻¹ of energy barrier, following the C oxidation path as shown in Fig. 2.

Fig. 2 Scheme of studied pathways for DRM on Rh(111).

On the other hand, CH* intermediates can also be oxidized by O* (CH* + O* → CHO*) before dissociated into C* + H*, overcoming an energy barrier of 156 kJ mol⁻¹ to produce HCO*, through the CH oxidation path (Fig. 2). CO* is subsequently obtained by HCO* dissociation into CO* + H*. The CH/C oxidation paths lead to the main DRM products, CO_(g), and H_2(g), with the CH oxidation pathway being slightly more favorable per barely 11 kJ mol⁻¹ of energy.

The H* produced by the methane activation/cracking can react with O* to form OH*, the precursor for the water formation, which can be alternatively obtained by the COOH* dissociation into CO* + OH*. The OH* formation involves practically the same energy barrier from the COOH* dissociation (119 kJ mol⁻¹) than from the reaction of O* + H* (120 kJ mol⁻¹), resulting in a lack of understanding about the most favorable pathway to form the DRM products by means only DFT calculations. Subsequently, OH* reacts with H* to produce water, a secondary product of the DRM, through the competitive RWGS path shown in Fig. 2. Note that the differences between different pathways studied here for the overall DRM mechanism are barely appreciable; thus, the conclusions about the main route operating over DRM on Rh(111) cannot be drawn from only DFT energies (see ESI† Fig. S4), and kinetic Monte Carlo simulations are needed. The energies presented in this section correspond to those used in the kinetic model, as described in section 2.2.

3.2. Kinetic model

kMC simulations are carried out under the typical operating reaction conditions,¹⁶i.e., pressures of 1 bar and temperatures from 700 K to 825 K, to obtain further insight into the mechanism of DRM on Rh(111). These simulations are performed for a CO₂/CH₄ mixture with a 1 : 1 ratio. First, we analyzed the event frequencies of the elementary steps and the turnover frequency (TOF), calculated as the ratio between the number of molecules of products (CO or H₂) over time and the number of sites, to estimate the catalytic activity and the dominant elementary steps controlling the reaction. Each simulation has been allowed to achieve a steady state in which no variations of the TOF are observed, except for small fluctuations resulting from the stochastic nature of the kMC method (ESI† Fig. S3). The simulation time necessary for all temperatures to reach the steady state is about 1 or 3 seconds.

The event frequencies of the elementary steps under different operating temperatures and at steady state are shown below.

As presented in Fig. 3, the CO₂ activation to CO* and O* can lead to two adsorption sites for CO, namely fcc and hcp sites, after its activation through transition-states with comparable energies. Nevertheless, the CO₂ activation leading to the CO* fcc has a higher occurrence frequency than that leading to the CO* hcp, possibly due to the competition for the hcp sites with C* and CH*, two of the most abundant species on the surface, converting the in the most feasible route to activate CO₂. Fig. 3 shows that the production of the DRM primary products on Rh(111) is dominated by the HCO* formation/activation through the CH oxidation pathway, according to what is suggested from the analysis of the DFT energy profiles. The frequency of CO* formation through the C oxidation pathway is about 3 orders of magnitude smaller than that through the CH oxidation pathway. Thereby, the CH oxidation path becomes the primary route controlling the formation of DRM reaction products. Still, the C oxidation pathway also contributes to forming CO, but to a minor extent. This role of both oxidation pathways in the DRM cannot be inferred from the analysis of energy profiles. The literature has also claimed the OH* molecules as key intermediates, but our results suggest that O* is the oxidant of the C/CH* species since OH* is not present on the surface.

Fig. 3 Event frequency per site for the dry reforming of methane reaction catalyzed by Rh(111) under different temperatures: a) T = 700 K; b) T = 825 K.

Regarding the RWGS pathway, the observed COOH* formation is practically negligible under any studied conditions (see ESI† Fig. S2 for all temperatures), as suggested by the DFT calculations, in which the direct CO₂ activation was clearly favored over the hydrogen-assisted one. Since only one reaction pathway, O* + H*, involves species produced from the direct activation of both reactants (O* from CO₂ and H* from CH₄), contributes to the RWGS reaction (O* + H* → OH*, the main precursor for forming water, see Fig. 3) water formation occurs but to a small extent, so it does not affect the final ideal H₂/CO ratio equal to 1. The water is formed only through the OH* hydrogenation (OH* + H* → H₂O*) since the disproportionation of two neighboring OH* species (OH* + OH* → H₂O* + O*) does not occur. This is probably due to the shallow coverage of OH* as the energy barrier for the OH* disproportionation is much lower than the OH* hydrogenation one: 44 vs. 78 kJ mol⁻¹, respectively.

Since the formation of species is relevant in designing DRM catalysts where coke formation is one of the biggest challenges to overcome, it is important to consider reactions between carbon atoms further than solely C/CH* oxidation reactions or even to produce CO*. However, the formation of CH₂O* via the oxidation path does not affect the kinetics and, therefore, is not included in our analysis (see ESI† Fig. S9). The carbon species are not directly comparable to coke formation on the surface. Still, they can be considered precursors in the carbon deposition process that could eventually lead to coke formation on the catalyst surface, resulting in deactivation. Therefore, we studied the formation of as a valuable indicator of the potential impact of carbon deposition over catalytic activity, and we found that the system is still active after including the C* migration and formation (see ESI† Fig. S8). The significant coverage of species found on the Rh(111) surface is not poisoning the surface; therefore, during the kMC simulations, there is no evidence of deactivation. The TOF (CO) decreased from 323 to 194 molecule per site per second, as expected when including a reaction that also consumes C*. However, the ratio between the main products is equal to 1. Therefore, we can conclude that the formation of carbonaceous species does not block the catalytic activity of Rh(111) for DRM.

The efficiency of the Rh(111) surface as a catalyst of the DRM reaction is described by the turnover frequency (TOF), defined here as the number of products formed (CO or H₂) per number of sites and time in seconds. For each simulation, the TOF is extracted as the slope of the number of CO (or H₂) molecules produced as a function of time divided by the number of sites after the steady state is reached (see ESI† Fig. S6).⁴⁶ Thus, different kMC simulations were run at different temperatures, revealing that increasing the temperature has no other effect than increasing the occurrence frequency for all events in a homogenous way, except for the frequency of the water formation step since it decreases slightly as the temperature increases. The overall TOF increases as the temperature increases, from 16 at 700 K to 65 at 825 K (see ESI† Table S2). These results suggest that the Rh(111) surface is less active than the Ru(0001) surface for the DRM reaction when comparing both systems using the same level of theory.²⁷

The increase in the catalytic activity (TOF) under steady-state conditions is not remarkable in the range of temperatures studied. In the same way, the coverage obtained at different temperatures does not present significant variations. The surface's main species are the C* and CH* intermediates. Afterward, the minor species on the surface are H* and from higher to lower abundance, respectively. Only the significant coverages, coverages ≥ 0.1, are plotted. As shown in Fig. 4, the most abundant species on the surface (C* and CH*) correspond to the most strongly adsorbed species on the Rh(111) surface, which fluctuate in time and compete for the occupation of the same hcp sites. At higher temperatures, the same trends are observed with a slight increase in C* coverage while CH* coverage decreases slightly. It is worth mentioning here the need to include the lateral interactions CH–CH*, C–C*, CH–C*, and C–O* since otherwise, the surface is saturated, causing a rapid decay of the catalytic activity. After including the lateral interactions between the most abundant species on the surface, it is observed that the catalytic activity is maintained for long simulation times. No evidence of deactivation was found during the kMC simulations until 10 seconds of total simulation time, and the TOF is already constant after 1–3 seconds. This contrasts with the faster convergence of the TOF in the Ru(0001) catalyst below 0.1 seconds, as found in our previous multiscale modeling study.²⁷ The main species covering the Ru(0001) surface were O*, C*, and CH*, with the latter having a lower coverage than those observed on the Rh(111) surface. Nevertheless, the fact that the coverage for both surfaces fluctuates over time probably allows the catalytic activity not to decay. However, the TOF achieved by Ru(0001) (2.57 × 10⁴ s⁻¹) was much higher than that achieved by Rh(111) at the steady-state conditions (3.23 × 10² s⁻¹) under reforming conditions (973 K, 1 bar; see ESI† Fig. S5 for Rh(111) kMC simulation at 973 K). This result agrees with an experimental work where Ru-based catalysts were more active than Rh-based ones.⁴⁷

Fig. 4 Surface coverage evolution, expressed as the fraction of active sites occupied per adsorbate, among the main species present on the surface during DRM reaction at 700 and 825 K.

3.2.1. Kinetic analysis of activity of DRM reaction. Several studies have evaluated a kinetic analysis of the CO₂ reforming of methane to provide fundamental kinetic parameters such as TOF, apparent activation energy, and partial pressure dependencies.⁴⁸ The reaction kinetics depends on the metal, as evidenced in the literature, but on Rh catalysts (Rh/SiO₂), experimental data suggest that CH₄ dissociation is the rate-limiting step.⁴⁹ In other studies, the selected model presents the surface reaction of both adsorbed reactants, CH₄ and CO₂,⁵ or even the surface reaction of carbon species (CH_x) and O,⁵⁰ as the rate-determining step. Overall, the DRM kinetic model continues to be a topic of debate.

Here, with the established kinetic model in the previous section, it is possible to unravel the kinetic behaviors of the activity of DRM reaction over the model Rh(111) surface. First, to identify the rate-determining step (RDS) controlling the overall activity of the DRM process taking place on the Rh(111) surface, the concept of degree of rate control (DRC) originally proposed by Campbell was applied.^51–53 The DRC can determine how each elementary step i affects the overall rate defined as:

where r is the overall rate to the product of interest (i.e., TOF), and the partial derivate is taken fixing all other rate constants, k_j, for all other steps j ≠ i and the equilibrium constant for step i, K_i. To keep K_i constant, the forward and reverse rate constants for step i, must be modified by equal factors so that their ratio remains constant. Above, X_DRC,i equals the relative increase in the overall reaction rate per relative increase in the rate constant for step i (differentially), which ranges from −1 to 1. The larger the numeric value of X_DRC,i for a given step i, the higher is the influence of its rate constant over the overall reaction rate r. The more positive value of X_DRC,i means a more rate-limiting step; conversely, the more negative value implies that the step is more rate-inhibiting.⁵⁴ In this work, all elementary steps were analyzed, except the quasi-equilibrated steps of CH₂/CH activation, which are expected to behave as the CH₃ activation, and the formation/activation of COOH* intermediate, since it does not occur, as shown in Fig. 3. We vary the rate constant of the elementary step of interest to find their respective DRC and propose a rate-limiting step of the DRM reaction catalyzed by Rh(111) based on DFT + kMC simulations. We modify only the barrier energy to keep the equilibrium K_i constant by ±0.05 eV, around the calculated value at the DFT level. Then, new kMC simulations are run each time, and the respective TOF in the steady state is calculated. Finally, the X_DRC,i is obtained, in this case, via two different numerical methods: finite difference (FD) approximation averaging over two simulation runs, and a polynomial method where we fit the changes of TOF to a polynomial over five simulation runs. A polynomial is used here to approximate the function around an origin value by truncating a Taylor series expansion at the fourth-degree variable to obtain the DRC value with a small error. Although the Taylor's theorem is the central concept in the FD method, it is an approximation that can introduce larger errors compared to a Taylor series taking more terms. The function ln(TOF) is fitted to a polynomial to minimize the error, where the linear term represents its first derivative with respect to ln(k_i). This linear term is crucial since, as shown in eqn (7), it corresponds to the DRC value.

As the FD method requires only two simulation runs, it was used for all elementary steps evaluated. In contrast, the polynomial method was considered only for the step with the most significant DRC. Considering the stochastic nature of kMC simulations, which can introduce a large amount of noise, the rate constants were changed by a relatively large magnitude, until 5, as proposed by Hess et al., to avoid this issue.⁵⁵ The X_DRC,i values are reported in Table 1.

Table 1 Campbell's degree of rate control of the key elementary steps on DRM reaction catalyzed by Rh(111) at 700 K of temperature

Elementary step	DRC
	FD	Polynomial
CO2(g) + 2* → CO* + O*	1.09	0.95
	−0.09	—
	0.00	—
C* + O* → CO* + *	0.04	—
CH* + O* → HCO* + *	0.04	—
HCO* + * → CO* + H*	0.01	—
	0.01	—
O* + H* → OH* + *	0.00	—
OH* + H* → H2O* + *	0.01	—

---

The X_DRC,i was investigated only at 700 K since no significant differences were found for other temperatures. The results at the steady state show that the CO₂ direct activation process is the only RDS at the studied conditions (700 K and 1 bar). Since the X_DRC,i for this step is positive, the overall rate of CO production will increase if the barrier of CO₂ activation decreases. The main DRM products, CO and H₂, have an ideal H₂/CO = 1 ratio, which is the main goal for industrial purposes. Thus, the formation rate of both products is critically dependent on the CO₂ dissociation. The fact that X_DRC,i for CO₂ activation is slightly greater than 1 (1.09), shows that the kMC simulations are accurate enough to reproduce accurate insights into the kinetics of the mechanism, considering the approximations of the numerical method used to solve eqn (7) and the stochastic nature of kMC simulations. The methane activation is the second step that contributes slightly to the degree of rate control but with a small inhibiting effect. For all other evaluated steps, the X_DRC,i is approximately zero.

To know more accurately the value of X_DRC,i, we employed the polynomial method for the rate-limiting step controlling the catalytic activity. For this method, 5 simulations were run only for the transition state of the CO₂ activation, varying the energy barrier from ±0.05 eV until ±0.15 eV; in other words, the rate constant variation was up to . Then, the changes of TOF were fit to a polynomial (ESI† Fig. S7), where the linear term is the desired derivative. The result of the polynomial method is shown in Table 1 in the second column. The result for the CO₂ activation step is close to that obtained through the finite difference approximation, as expected due to the small obtained coefficients in the polynomial beyond the first term, revealing that the CO₂ activation is the rate-limiting step, in contrast with the experimental finding over Rh/SiO₂.⁴⁹ In fact, at the earliest stages during all kMC simulations, we observed that the H₂ product is much more abundant than CO, suggesting that CH₄ activation occurs faster on the surface than CO₂ activation. Therefore, the CH_x species become available on the surface with a relative abundance higher than those from CO₂ activation in a shorter time. Thus, the relevant reactions between CH₄ cracking (CH_x) products and O* to produce the second stoichiometry CO molecule and the key removal of adsorbed C* on the surface to continue the reaction depend strictly on CO₂ activation. This DRC analysis over CO₂ activation was also carried out following the changes in the H₂ production and the total TOF, obtaining practically the same value of DRC (see ESI† Fig. S7). Thus, we only report the TOF based on CO for all cases here.

3.2.2. Apparent activation energy. Another key kinetic parameter of the DRM reaction is the apparent activation energy, which ranges from 29 to 306 kJ mol⁻¹, depending on the catalyst used.⁴⁸ Here, the apparent activation energy of the overall reaction was computed by performing different kMC simulations at six different temperatures, namely 700, 720, 740, 775, 800, and 825 K, according to the experimental conditions of the DRM reaction when catalyzed by Rh-based catalysts.¹⁶ Plotting the steady-state activity data as an Arrhenius plot, i.e., as the logarithm of TOF versus 1000/T, a straight line was obtained whose slope corresponds to the apparent activation energy. Fig. 5 shows an Arrhenius plot of the rates of CO formation of DRM reaction in a CO₂/CH₄ 1 : 1 reaction mixture with a pressure total of 1 bar. The calculated activation energy associated with the plot in the figure is 53 kJ mol⁻¹. Unfortunately, to our knowledge, the experimental apparent activation energy of the DRM catalyzed by the Rh(111) has not been reported so far. There is, however, an experimental study where various Ni-based bimetallic catalysts for DRM are compared at low temperatures, resulting in the Ni–Rh/Al₂O₃ catalyst being the most active one exhibiting the highest conversions of CH₄ and CO₂ due to a low energy activation of 54 kJ mol⁻¹. This activation energy was calculated from the Arrhenius plot at 1 bar of pressure and temperature ranging from 673 to 873 K, and it is very close to our apparent activation energy calculated here for Rh(111): 53 kJ mol⁻¹.

Fig. 5 Arrhenius plot of the DRM TOF (in s⁻¹ per site) in the temperature range of 700–825 K and at 1 bar of pressure.

The previous section showed that the DRC criteria provide useful information about the system's chemistry and can be employed to obtain guidance about the parameters that critically determine the overall catalytic activity. Here, we use this criterion to explain the apparent activation energy of the network (E_app), as was demonstrated in the literature, using the following approximation:⁵⁶

where the apparent activation energy (E_app) is approximately given by an “average” of the activation energies of all elementary processes (ΔE^≠_i) weighted with their DRC criteria X_DRC,i. Thus, if there is only one process with a much higher DRC than the rest, the E_app of the whole reaction only depends on the activation energy of that step. Such a situation is not general, as multiple processes could control the overall rate, but this seems true for the DRM catalyzed by Rh(111). Note that the apparent activation energy calculated from the Arrhenius plot is close to the energy barrier of the CO₂ direct dissociation step (i.e., 47 kJ mol⁻¹), thus confirming that this is the rate-determining step, as expected from the DRC analysis and the E_app approximation.

4. Conclusions

We explored the kinetics of dry reforming of methane catalyzed by Rh(111) by combining DFT calculations with kinetic Monte Carlo simulations in the temperature range of 700–825 K and 1 bar of pressure. Our study encompassed 38 elementary reactions using Langmuir–Hinshelwood and Eley–Rideal formalisms, including the competitive RWGS reaction. The kinetic model considered CH₄ and CO₂ dissociative adsorption, surface chemical reactions between intermediates, CO_2(g) reactions with adsorbed H*, and gas-phase product desorption. Our findings revealed that the CH oxidation pathway is preferred to produce CO_(g), surpassing the C oxidation pathway. O* played a crucial role in oxidizing C/CH* species. At steady-state conditions, abundant C* and CH* were found, and lower concentrations of H*, , and O*, indicating available O* for oxidizing carbon species even after prolonged simulation. Including lateral interactions is essential to prevent surface saturation with C*, CH*, and O*, affecting reaction kinetics and barrier energies.²⁰ Our kMC simulations confirmed the dominance of the RWGS reaction via CO₂ direct dissociation, attributed to its lower activation barrier compared to COOH* formation. In any case, the water formation has a negligible effect on the final H₂/CO ratio.

We employed the DRC criteria to identify the rate-determining step, revealing that CO₂ activation is the only rate-limiting step. Therefore, the calculated E_app approximated the energy barrier for CO₂ dissociation. Our work underscores the significance of multiscale modeling in understanding complex catalytic systems, providing insights into turnover frequency, product distribution, main reaction pathways, surface coverage, and crucial atomistic steps. This approach enhances our understanding of catalytic processes, especially in heterogeneous systems, overcoming the limitations of Gibbs energy profiles. Moreover, kMC modeling can account for multiple adsorbed species, lateral interactions, and distinct active sites, allowing for more accurate modeling than via mean-field microkinetic modeling since the latter simplifies the surface species' dynamics and behavior.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This research was funded by the Spanish “Ministerio de Ciencia e Innovación” (PGC2018-100818-A-I00 and PID 2021-128416NB-I00). The authors thank the undergraduate students Laura Fernández and Trinité Perant for their contributions to the computational study of the absorption of the DRM key intermediates and the CH₄ activation step.

References

1. C. Le Quéré, J. I. Korsbakken, C. Wilson, J. Tosun, R. Andrew, R. J. Andres, J. G. Canadell, A. Jordan, G. P. Peters and D. P. van Vuuren, Nat. Clim. Change, 2019, 9, 213–217.
2. P. Friedlingstein, M. O. Sullivan, M. W. Jones, R. M. Andrew, L. Gregor, J. Hauck, C. Le Quéré, I. T. Luijkx, A. Olsen, G. P. Peters and W. Peters, Earth Syst. Sci. Data, 2022, 14, 4811–4900.
3. D. Chery, V. Lair and M. Cassir, Front. Energy Res., 2015, 3, 1–10.
4. D. Pakhare and J. Spivey, Chem. Soc. Rev., 2014, 43, 7813–7837.
5. M. M. B. Quiroga and A. E. C. Luna, Ind. Eng. Chem. Res., 2007, 46, 5265–5270.
6. K. Aasberg-Petersen, J. H. Bak Hansen, T. S. Christensen, I. Dybkjaer, P. S. Christensen, C. Stub Nielsen, S. E. L. Winter Madsen and J. R. Rostrup-Nielsen, Appl. Catal., A, 2001, 221, 379–387.
7. M. Akri, S. Pronier, T. Chafik, O. Achak, P. Granger, P. Simon, M. Trentesaux and C. Batiot-Dupeyrat, Appl. Catal., B, 2017, 205, 519–531.
8. Q. Zhang, M. Akri, Y. Yang and B. Qiao, Cell Rep. Phys. Sci., 2023, 4, 101310.
9. G. S. Gallego, C. Batiot-Dupeyrat, J. Barrault, E. Florez and F. Mondragón, Appl. Catal., A, 2008, 334, 251–258.
10. O. Muraza and A. Galadima, Int. J. Energy Res., 2015, 39, 1196–1216.
11. J. R. Rostrup-Nielsen and J. H. Bak Hansen, J. Catal., 1993, 144, 38–49.
12. L. Foppa, T. Margossian, S. M. Kim, C. Müller, C. Copéret, K. Larmier and A. Comas-Vives, J. Am. Chem. Soc., 2017, 139, 17128–17139.
13. A. Tarasov, H. Düdder, K. Mette, S. Kühl, K. Kähler, R. Schlögl, M. Muhler and M. Behrens, Chem. Ing. Tech., 2014, 86, 1916–1924.
14. Z. Liu, D. C. Grinter, P. G. Lustemberg, T. D. Nguyen-Phan, Y. Zhou, S. Luo, I. Waluyo, E. J. Crumlin, D. J. Stacchiola, J. Zhou, J. Carrasco, H. F. Busnengo, M. V. Ganduglia-Pirovano, S. D. Senanayake and J. A. Rodriguez, Angew. Chem., Int. Ed., 2016, 55, 7455–7459.
15. M. Usman, W. M. A. Wan Daud and H. F. Abbas, Renewable Sustainable Energy Rev., 2015, 45, 710–744.
16. J. A. Lee, Y. Bae, K. Hong and J. Hong, Int. J. Energy Res., 2022, 46, 11228–11249.
17. A. Bruix, J. T. Margraf, M. Andersen and K. Reuter, Nat. Catal., 2019, 2, 659–670.
18. N. López, N. Almora-Barrios, G. Carchini, P. Błoński, L. Bellarosa, R. García-Muelas, G. Novell-Leruth and M. García-Mota, Catal. Sci. Technol., 2012, 2, 2405–2417.
19. M. Stamatakis, J. Phys.: Condens. Matter, 2015, 27, 1–28.
20. M. Pineda and M. Stamatakis, J. Chem. Phys., 2022, 156, 1–30.
21. 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, p. 315.
22. L. V. Lutsevich, O. A. Tkachenko and V. P. Zhdanov, Langmuir, 1992, 8, 1757–1761.
23. J. Nielsen, M. D'Avezac, J. Hetherington and M. Stamatakis, J. Chem. Phys., 2013, 139, 1–13.
24. M. Stamatakis and D. G. Vlachos, J. Chem. Phys., 2011, 134, 1–13.
25. H. Prats, S. Posada-Pérez, J. A. Rodriguez, R. Sayós and F. Illas, ACS Catal., 2019, 9, 9117–9126.
26. X. Li and L. C. Grabow, Catal. Today, 2022, 387, 150–158.
27. E. Díaz López and A. Comas-Vives, Catal. Sci. Technol., 2022, 12, 4350–4364.
28. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
29. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
30. Z. Zhang, S. S. Wang, R. Song, T. Cao, L. Luo, X. Chen, Y. Gao, J. Lu, W. X. Li and W. Huang, Nat. Commun., 2017, 8, 1–10.
31. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
32. J. Wellendorff, K. T. Lundgaard, A. Møgelhøj, V. Petzold, D. D. Landis, J. K. Nørskov, T. Bligaard and K. W. Jacobsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 32–34.
33. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
34. C. Papadopoulou, H. Matralis and X. Verykios, Catalysis for Alternative Energy Generation, Springer N, 2012.
35. L. Foppa, M. C. Silaghi, K. Larmier and A. Comas-Vives, J. Catal., 2016, 343, 196–207.
36. D. A. McQuarrie, Statistical mechanics/Donald A. McQuarrie, 1975.
37. M. Stamatakis and D. G. Vlachos, ACS Catal., 2012, 2, 2648–2663.
38. M. Stamatakis, Zacros Software, https://zacros.org, 2013.
39. K. A. Fichthorn and W. H. Weinberg, J. Chem. Phys., 1991, 95, 1090–1096.
40. K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 1–17.
41. K. J. Laidler, Chemical Kinetics, Harper & Row, 1987.
42. M. T. Darby, S. Piccinin and M. Stamatakis, in Physics of Surface, Interface and Cluster Catalysis, 2016, pp. 4.1–4.38.
43. R. Prins, Top. Catal., 2018, 61, 714–721.
44. E. K. Rideal, Math. Proc. Cambridge Philos. Soc., 1939, 35, 130–132.
45. W. H. Weinberg, Surf. Sci. Lett., 1983, 128, 12–14.
46. D. W. Goodman, D. E. Peebles and J. M. White, Surf. Sci. Lett., 1984, 140, 239–243.
47. A. Androulakis, I. V. Yentekakis and P. Panagiotopoulou, Int. J. Hydrogen Energy, 2023, 48, 33886–33902.
48. M. C. J. Bradford and M. A. Vannice, Catal. Rev.: Sci. Eng., 1999, 41, 1–42.
49. H. Y. Wang and C. T. Au, Appl. Catal., A, 1997, 155, 239–252.
50. T. Osaki, H. Masuda and T. Mori, Catal. Lett., 1994, 29, 33–37.
51. C. T. Campbell, Top. Catal., 1994, 1, 353–366.
52. C. T. Campbell, J. Catal., 2001, 204, 520–524.
53. C. T. Campbell, ACS Catal., 2017, 7, 2770–2779.
54. W. Xie, J. Xu, Y. Ding and P. Hu, ACS Catal., 2021, 11, 4094–4106.
55. F. Hess and H. Over, ACS Catal., 2017, 7, 128–138.
56. H. Meskine, S. Matera, M. Scheffler, K. Reuter and H. Metiu, Surf. Sci., 2009, 603, 1724–1730.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cy01546g

---