Unraveling the prominent role of the Rh/ZrO2-interface in the water-gas shift reaction via a first principles microkinetic study

The industrially important water–gas shift (WGS) reaction is a complex network of competing elementary reactions, where the catalyst is a multicomponent system consisting of distinct domains. Herein we have combined density functional theory calculations with microkinetic modeling to explore the active phase, kinetics, and reaction mechanism of the WGS over Rh/ZrO2. We have explicitly considered the support, metal, and their interface and find that the Rh/ZrO2 interface is far more active towards WGS than Rh(111) facets, which are susceptible to CO poisoning. CO2 forming on the zirconia support rapidly transforms into formate. These findings demonstrate the central role of the interface in the water–gas shift reaction, and the importance of modeling both the support and the metal in bifunctional systems.

simple five step Langmuir-Hinshelwood kinetic model has been previously used to represent experimental data on the WGS over a Rh/ZrO 2 catalyst. 27 The model does not make assumptions about the possible intermediates of the surface reaction steps, although hydroxyl was suggested to be one of the important reaction intermediates in an earlier study. 22 In that study, the Mars-van-Krevelen type redox mechanism was ruled out based on the observation that rhodium was in a metallic form under WGS reaction conditions. Hydroxyl, carbonates, and carboxylate were suggested to be possible reaction intermediates, while formates were observed on the zirconia support of the catalyst at temperatures below 623 K 22 .
In general, multiple reaction mechanisms have been put forward for the water-gas shift reaction on oxide supported metal catalysts, and the dominant mechanism depends on both the nature of the catalyst and the support. It is commonly suggested that WGS prefers the so-called "redox" mechanism over reducible oxide supports. 1,28,29 In this Mars-van-Krevelentype mechanism the oxide support provides lattice oxygen for the oxidation of CO into CO 2 leaving an oxygen vacancy behind. The oxide is subsequently reoxidized by water which dissociates into the oxygen vacancy producing hydrogen in the process. Irreducible oxides, such as zirconia, 30 are thought to be incapable of providing lattice oxygens for the redox reaction, and therefore alternative "associative" mechanisms have been put forward for catalysts supported on those materials, 12,15 as well as for CeO 2 supported systems. 10,11 In associative pathways, CO reacts on the catalyst surface with dissociated water to produce an intermediate (e.g. carboxyl or formate) species which then dissociates into a hydrogen atom and CO 2 . The nature of this intermediate has been under debate. 10,11,15,[31][32][33][34][35][36][37][38] Several computational studies on different unsupported metal surfaces suggest that carboxyl is the reaction intermediate for the associative mechanism. [31][32][33][34] However, many experimental groups have observed formate on the oxide support under WGS reaction conditions and have argued it to be the main reaction intermediate. 10,11,35,37,38 Formate is said to be formed from the reaction between CO and a "terminal" OH group bound oxygen down to a single metal cation of the oxide, generating a bidentately bound formate. 11,15,36,38,39 The reaction is unlikely to occur 3 in a single low barrier step due to the many bonds that are required to be broken and formed simultaneously. In contrast, carboxyl can be formed in a single step from CO and OH.
Another pathway that also does not require the formation of oxygen vacancies on the support is the "surface redox" mechanism, in which water fully dissociates to form H 2 and atomic oxygen, which can then oxidize CO into CO 2 . 40 Complete water dissociation can take place via two different pathways. After the initial partial dissociation of water into OH and H, atomic oxygen can be formed from direct dissociation of OH, or from disproportionation reaction between two OH groups to form H 2 O and O. As the surface redox pathway does not require the reduction and oxidation of the oxide support, it can, in principle, operate on the metal alone, and has been included in many DFT studies on the WGS on unsupported metal surfaces. [31][32][33][41][42][43] Surface redox was originally thought to be the main reaction mechanism on the Cu(111) surface, 40 however it was later shown in a DFT study that WGS proceeds through a carboxyl mechanism on Cu(111) instead. 31 The surface redox route has been recently shown to be the dominant mechanism for Ni/Al 2 O 3 in a combined computational and experimental study. 18 Computational studies typically employ slab models to perform calculations for the WGS over metal catalysts, [31][32][33][41][42][43] however, since the slab contains only metal sites, the models completely ignore the possible role of the underlying oxide supports, and by extension, the catalyst-support interface. The forward and reverse WGS reactions over rhodium have also been previously studied employing slab models of the Rh(111)-facet. 33,42,44 These studies employed only DFT calculations, and no solid understanding was reached about the dominating mechanism for the WGS on Rh(111). The role of formate was only explicitly considered in Ref. [ 33], where it was stated that the formate pathway was energetically competitive with the redox and carboxyl pathways. That study suggests that formate is the most abundant experimentally observed intermediate because its formation barrier is lower than its dissociation barrier. Overall, the rationalization of obtained DFT results on WGS calls for further analysis, such as a microkinetic modeling, which allows resolving the dominant pathway 4 under varying reaction conditions. Furthermore, a catalyst model including both the metal and the support might be needed to better represent the material of interest.
While catalyst-support interfaces are in principle straightforward to model, constructing a computationally feasible model including necessary details is difficult. In cluster models, the metal catalyst consisting of up to several tens of atoms is placed on an oxide slab. 16,18,[45][46][47][48][49] This can be a useful technique for modeling metal clusters and small particles, but is unsuitable for modeling nanoparticles in the range larger than e.g. 3 nm, active for WGS on Rh/ZrO 2 . 22 Recently, computational studies have employed oxide supported metal nanorods to model the metal-oxide interface. 8,9,17,19,[50][51][52] In these models, a 2 − 3 layers thick metal nanorod is placed on an oxide surface slab. The calculations are periodic along the length of the rod, meaning the rod is infinitely long. This construction can be taken to represent an actual nanowire with a cross section of a few atoms, or the perimeter of an arbitrarily large nanoparticle.
First principles microkinetic modeling is often used to unravel the complex reaction network and relate the zero Kelvin and zero pressure DFT results to the real catalytic process taking place at reaction conditions. 4,53,54 For WGS reaction, there exists DFT based microkinetic models for unsupported metal surfaces, 31,32 but as explained above, those calculations exclude the potentially important contribution from the support. Recently, DFT based WGS microkinetic models have been constructed for oxide supported metal structures, namely Ni/Al 2 O 3 , 18 Au/TiO 2 , 19 , Cu/ZrO 2 , 51 and Au/MgO. 8 All of these studies highlight the metal-oxide interface to be of significant importance to the WGS activity for the materials studied, information which would not be attainable without an explicit interface model.
For the Ni 55 /γ-Al 2 O 3 cluster model, the surface redox reaction was shown to occur at the metal-support interface sites, with OH dissociation as the rate limiting step. 18 For the gold catalyst supported on the irreducible MgO, represented by a nanorod model, WGS was shown to proceed via the carboxyl mechanism. 8 The carboxyl formation step was found to be rate limiting. In contrast, when gold is supported on reducible TiO 2 , WGS proceeds 5 through the Mars-van-Krevelen-type redox reaction with water dissociation as a rate limiting step, as demonstrated using a nanorod model for Au/TiO 2 (110). 19 Comparing the MgO and TiO 2 supported gold systems reveals how great an impact the choice of support has on the mechanism. A Cu/m-ZrO 2 (212) rod model was used to study the reactivity of copper nanowires supported on zirconia. 51 It was found that although the interfacial sites do not directly take part in the WGS mechanism, the ZrO 2 support facilitates the WGS by modifying the Cu sites in the immediate vicinity of the oxygen-rich interface. The reaction was found to occur through the carboxyl mechanism, with water dissociation as the rate limiting step. All of these studies highlight the fact that depending on the identity of the metal and support, the WGS can occur through vastly different reaction mechanisms, the promoting effect of the interface being complex and varied.
In this work, we undertake a multiscale computational study of WGS over m-zirconia supported rhodium nanoparticles to elucidate the reaction mechanism and to understand how different constituents of a catalyst, the support, "bulk" nanoparticle, and the catalyst/supportinterface contribute to overall activity and what is the active phase of a catalyst. We use three model systems to capture the bifunctional nature of the Rh/ZrO 2 catalyst and include a complete set of elementary reaction steps to construct a WGS reaction network.
We begin by presenting DFT results for the thermodynamic and kinetic properties of the reaction. With first principles microkinetic modeling, we analyze the reaction mechanism and the contribution of the components of the catalyst towards CO conversion. We close by comparing our findings to existing experimental observations and highlight the complexity of water-gas shift chemistry, the challenges in microkinetic modeling, and the importance of different sites and parts of a catalyst to form overall understanding of the reaction at hand. 6

Computational Details
All DFT calculations reported herein were performed using the GPAW 1.1.0 software. [55][56][57][58] The exchange and correlation effects were described within the generalized gradient approximation employing the Perdew-Burke-Ernzerhof (PBE) functional 59,60 and the core electrons of elements were treated with the PAW 61 setups with the frozen core approximation. The Brillouin zones for the rhodium and zirconia surfaces were sampled using a (6 × 6 × 1) and a (4 × 4 × 1) Monkhorst-Pack mesh of k-points, respectively, while the interface and the gas phase species were treated at the Γ-point.
The transition states (TS) were located using a climbing image automated nudged elastic band (CI-AutoNEB) method [62][63][64][65] or constraint optimization, and they were confirmed by the presence of a single imaginary vibrational mode along the reaction coordinate. For full list of imaginary frequencies, please refer to the supporting info (SI). The vibrational analysis was performed for the reaction intermediates on all three model systems in order to calculate zero-point energies (ZPE). In each case, the analysis was performed only for the adsorbate atoms, and frequencies less than 100 cm −1 in energy were omitted from the ZPE calculations. Reaction rates for surface reactions were computed using harmonic transition state theory. 66 For non-activated adsorption steps, the rates were calculated using kinetic gas theory (see SI). Entropy of adsorbed species was calculated using the harmonic approximation as implemented in ASE, 57 while gas-phase entropies were taken from the NIST database. 67 We thoroughly screened possible adsorbate structures at the interface with the LCAO 68 DZP basis allowing them to relax until the maximum residual force was below 0.1 eV·Å −1 .
The most promising interface structures as well as rhodium and zirconia surfaces were then optimized in the grid basis with a 0.2Å maximum grid spacing. In grid basis, the structures were allowed to relax until the maximum residual force was below 0.05 eV·Å −1 .
The interface was modeled with a 40-atom Rh nanorod covalently bound to a 2ML-thick 7  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 m-ZrO 2 (111) surface with the Rh(111) facet of the front edge of the rod facing towards zirconia to mimic the perimeter of a large nanoparticle, please see Fig. 1c. The m-ZrO 2 (111) slab was chosen to represent the support since monoclinic zirconia is the most stable polymorph at typical WGS conditions, 69 and (111) the most stable facet. [69][70][71] The size of the unit cell of the ZrO 2 slab was (2x2), which was chosen to minimize lattice mismatch between the oxide and the Rh-rod. For the size of the chosen unit cell, the lattice mismatch between bulk rhodium and m-ZrO 2 (111) is -1.0 % and 6.8 % along the axes a and b (see Figure   1c), respectively. Since the calculations have to be periodic along the length of the rod, the alignment where the lattice mismatch is smallest, -1.0 %, was chosen. This means that in the  52 We also observe a greater ligand effect as compared to the Rh/MgO-interface, making the overall change in the d-band center more negative in our case. 52 To eliminate the interaction between an adsorbate and the back side of the next periodic image of the rod, adsorbates were only placed on oxide sites close to the Rh/ZrO 2 interface.
On the rod, adsorbates were placed on the first row of Rh atoms. The most stable extended rhodium (111) surface (Fig.1a) was modeled with a (2 × 2) slab with a lattice constant of 3.857Å. The slab was four atomic layers thick and the two bottom most layers were kept 8

Microkinetic Model
In a microkinetic model the pressure and surface coverage of each species as a function of time is obtained by solving the reaction rate equations. Rate constants are based on DFT-computed activation free energies and pre-factors are taken to be equal to k b T h . The complete reaction network and details of the microkinetic model are given in the SI. For brief, we used an in-house code to solve the reaction rate equations for a continuous stirred tank reactor. We solve the dynamical equations explicitly as a function of time because this gives important information on the transient time dependent dynamics of the complex system. Furthermore, this approach provides information on the reaction mechanism and how different species react, form, and disappear with time. The simulation time was set to 10000 s for all runs to ensure that the system had enough time to reach dynamic steady state. Catalyst mass, BET, Rh loading, and particle diameter were set at 100 mg, 100 cm 2 /g, 0.5 wt.%, and 5 nm, respectively, to mimic experimental values. For analyzing the results, we performed rate control analyses (for more details see SI). 72 From the time-dependent fluxes for (ṅ i (t)) reactants and products, we calculated the conversion as where in and out refer to in-put and out-put feeds while P is the pressure. We approximate the full conversion to be the time-independent steady-state (SS) conversion X SS . Reaction order with respect to the reactant and product gases was also determined by varying the partial pressure of one gas while keeping the others constant. For full details please see SI.

Results
We start with by presenting DFT results for a WGS reaction thermodynamics and kinetics over three different components of a Rh/ZrO 2 catalyst. Then we analyze and compare the DFT-computed reaction mechanisms collectively. Finally, we report predicted CO conversions and surface coverages from the microkinetic analysis.

Results from DFT calculations
We have considered the three most commonly proposed reaction pathways for the WGS: Associative carboxyl route, associative formate route, and the redox route.

Water Dissociation and Disproportionation
Regardless of the overall reaction mechanism, the WGS reaction always begins with at least partial water dissociation leading to the formation of hydrogen and a hydroxyl species followed by either a reaction with CO or a complete water dissociation to atomic hydrogen and oxygen. Figure 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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58 59 60  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 Complete water dissociation is a possible source of atomic oxygen for the surface redox reaction pathway. The reaction is endothermic on both m-ZrO 2 (111) and at the interface, with an activation energy of 1.34 and 1.00 eV, respectively. On Rh(111), the reaction is exothermic by −0.28 eV and has an activation energy of 0.85 eV. Alternatively, O atoms can also be generated on a catalyst surface by the disproportionation reaction between two OH species, where one OH species donates a proton to the other, forming atomic oxygen and water. On the Rh(111) surface, the reaction is mildly exothermic by −0.19 eV and occurs spontaneously. It has been previously shown that the disproportionation reaction proceeds spontaneously also on the extended Pt(111) surface and at the Ni/Al 2 O 3 -interface. 18,32 On zirconia, no stable configuration for coadsorbed atomic oxygen and molecular water could be found, but atomic structure relaxation lead to the spontaneous formation of two adsorbed OH species instead, which indicates that disproportination is unlikely on bare ZrO 2 . At the interface, molecular water and atomic oxygen are stable when placed far away from each other, but spontaneously form two OH groups when brought close together, making disproportination unfeasible also at the interface. is slightly less exothermic than reported in literature. 39 As discussed above, atomic oxygen may be obtained from disproportionation or complete water dissociation. Its adsorption energy with respect to gas-phase water is +0.19 eV and +0.31 eV for the interface and rhodium systems, respectively. Prior to CO oxidation, coad-  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 most favorable on the zirconia, since the reaction takes place without a kinetic barrier and is exothermic by −0.25 eV. However, the difficulty of forming atomic oxygen on zirconia due to the high barrier of full water dissociation and OH disproportionation renders CO oxidation highly unlikely. The TS geometries on Rh(111) and the interface are both "L-shaped" with C-O bond lengths of 1.81Å and 1.80Å, respectively. The TS at the interface resembles the final state of the reaction, while on Rh(111) the TS is more similar to the initial structure.

CO adsorption and oxidation
Besides the difference in activation energies, the strength of CO 2 adsorption varies on the three surfaces; the molecule adsorbs exothermally at the interface and on pure zirconia by −0.68eV and −0.57eV, respectively, but is only weakly physisorbed on the Rh(111) by −0.01eV. This implies that once formed, the CO 2 molecule is likely to desorb immediately from the rhodium surface to the gas-phase, whereas on the zirconia and at the interface, it is available to react further with some other surface species like hydrogen.

Carboxyl formation and dissociation
Carboxyl is formed from coadsorbed OH and CO, and can have either a cis or trans geometry, the latter isomer being the more stable on all three model systems. on the Rh(111), while OH prefers the Zr cation. On pristine Rh and ZrO 2 surfaces, either CO or OH has to adsorb at an unfavorable site. At the interface, both metal and oxide sites are available to accommodate the adsorbates; CO is bound to a Rh atom and OH is bound on a Rh-Zr dual-site, leading to a more stable coadsorbed structure and higher activation energy than for the two pristine surfaces. The carboxyl dissociation step is most feasible at the interface, as it has a barrier of only 0.37 eV, and is exothermic by −0.28 eV. On zirconia and rhodium the reaction is endothermic by +0.76 eV and +0. 23 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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 60

Formate formation and dissociation
The discussion of formate formation and dissociation is divided into two parts. First, we compare the pathway on the rhodium and interface, which is followed by the discussion concerning formate formation on bare zirconia. As noted previously by other DFT studies on the WGS, 31,32,51 formate formation in a single step is improbable due to the adsorption geometries of formate, as well as its precursors OH and CO. Instead, CO can react with a hydrogen atom to form a formyl intermediate which can then be oxidized by atomic oxygen into formate. Figure 5 shows the PES of formate formation and dissociation on the Rh (111) surface and at the interface. CO hydrogenation to form a formyl species is endothermic, the reaction energy being +1.17 eV and +0.40 eV on rhodium and the interface, respectively.  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57   This is due to differences both in the transition state and final state geometries. In our case, formate dissociation is easier than its formation, indicating that the clean metal surface is unable to accumulate formate.
We note that the the reaction between CO 2 and H to form formate has a barrier of 0.68 eV and is exothermic by −0.63 eV at the Rh/ZrO 2 -interface. Interestingly, very similar results have been found for the Cu 38 /m-ZrO 2 (111)-interface, where formate acts as an intermediate for the CO 2 hydrogenation into methanol. 48 Since we find formate formation to be easier from WGS products than from reactants, it indicates that formate could potentially act as a reaction inhibitor.
Next, we address formate formation on the zirconia surface, where we were unable to locate a good transition state with a competitive activation energy for the formation of the formyl intermediate. The following two pathways were considered. In formyl formation from adsorbed CO and a hydrogen bound to a lattice oxygen, the CO inserts itself into the O-  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 H bond and forms a "lattice formyl". The transformation of "lattice formyl" to bidentate formates requires extraction of oxygen from the lattice. This is typically costly, since the reduction of zirconia is strongly endothermic if no rhodium atoms or tiny clusters are in the vicinity. 76 In the other reaction pathway, CO reacts with a terminal OH group, again inserting itself into the O-H bond, which leads to a reaction pathway for which the location of transition state fails. One of us has previously investigated a similar reaction pathway, and found that the direct reaction of CO with a terminal OH has an activation energy of at least 1.60 eV. 39 Thus, we suggest the following alternative route for formate formation from WGS products on bare zirconia, as shown in Figure 6. First, adsorbed CO reacts with a terminal OH to form a carboxyl intermediate, which can then dissociate into adsorbed CO 2 and atomic hydrogen as discussed above. The newly formed CO 2 can either desorb or further react with hydrogen to formate, which is strongly exothermic by −2.45 eV, and has an activation energy of only 0.02 eV. Since CO 2 desorption is endothermic by +0.57 eV, formate formation is preferred. Due to the strong exothermicity of the reaction, formate acts as a thermodynamic "sink", therefore, eventhough CO 2 can also react with hydrogen back to carboxyl and even further back to CO and OH, the equilibrium is shifted towards formate.  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 Comparison of mechanistic aspects of studied reaction pathways Based on our DFT results show that water activation is more favorable both on the zirconia support, and at the Rh/ZrO 2 -interface, than on Rh(111). This is in line with the previous experimental findings which suggest that the role of the oxide support is to activate water. [11][12][13]28 Furthermore, the Rh/ZrO 2 interface can provide OH and H species for the subsequent steps of the WGS, thus no diffusion from the bulk oxide to the interface is required.
Interstingly, the identity of the metal seems to play a role, since for the Cu/ZrO 2 -interface, H 2 O dissociation was previously found to be "too facile". 51 This leads to strong adsorption of OH species at the interface and high CO oxidation barriers, rendering the Cu/ZrO 2 -interface unreactive. 51 For the Rh/ZrO 2 -interface, all barriers in the carboxyl route are lower than those in the surface redox route, which is limited by the high barrier of full dissociation of water. In particular, the carboxyl formation step which consumes OH, is favored over OH dissociation. Formation of WGS products via the formate intermediate is also unlikely, since it also requires full water dissociation. Furthermore, once formed, formate is very stable, and decomposition to either HCO + O or CO 2 + H requires overcoming a barrier of 1.63 eV or 1.31 eV, respectively. This could lead to the blocking of active sites by formate. As formate can also be formed from CO 2 on the surface, this is a possible source for CO 2 inhibition of the WGS. This has been observed experimentally for the Rh/ZrO 2 catalyst, 22,27 and has been computationally suggested for other systems. 18,31 On the Rh(111) surface, the highest barrier in the carboxyl pathway is the carboxyl dissociation step, while for the redox reaction the highest barrier is CO oxidation. The activation energy of the carboxyl dissociation step is 0.30 eV lower than that of the CO oxidation step, which means that the associative carboxyl pathway is energetically less costly of the two. Availability of OH groups for carboxyl formation is an important factor for the reaction pathway, as the disproportionation reaction readily consumes OH to form oxygen  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 barrier of 1.44 eV, and the equilibrium is towards the backward reaction into CO + H. The rhodium surface is not expected to accumulate formate, since the dissociation barrier is lower than the formation barrier. Because CO 2 only weakly interacts with the metal surface, readsorption of CO 2 , and its reaction with hydrogen to formate is unlikely. In addition, the strong binding energy of CO on Rh(111) indicates that the surface will have a high CO coverage, which might lead to self-poisoning and preventing water adsorption.
The zirconia support is not expected to be active for the WGS, 15 although the carboxyl route is shown to be a plausible way to generate the products on the m-ZrO 2 (111) surface. One possible reason for the inactivity is that the formation of formate from CO 2 and hydrogen has a very low barrier of only 0.02 eV, and the reaction is strongly exothermic, essentially trapping the products on the surface. The complexity of the network prevents further analysis of the dominant pathway and comparison between the activity of different catalyst domains based on DFT results.

Microkinetic Analysis
To fully elucidate how the pathways discussed above compete under reaction conditions, a comprehensive microkinetic analysis is performed. Microkinetic-reactor model simulations predict consumption of reactants and formation of products as a function of time under given reaction conditions. Herein, we report the results from simulations over three different model systems: a flat Rh(111) surface, a Rh(111)/ZrO 2 (111)-interface, and an ideal ZrO 2 (111) surface. The number of sites in the individual models were estimated to mimic the ratios of the metal, support, and interface sites present in a real catalyst (see SI). For our chosen particle size of 5 nm, the ratio of interface to metal facet sites is 1:12. In the microkinetic analyses temperature was varied from 300 K to 800 K at 100 K intervals and total pressure was set at 2 bar but water to CO ratio was altered. The initial coverage of all intermediates was set to zero for all simulations. No a priori assumptions were made about the rate determining step or steady state gas composition and surface coverages. The model 21  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 remains thermodynamically consistent through detailed balance. The equilibrium ratio Kp Keq was evaluated for all simulations at the steady state conditions, and was found to be less than 0.1 in all cases, meaning that the steady state was far from the thermodynamic equilibrium for all simulations. For further details concerning the microkinetic model, please refer to the SI.
WGS on ZrO 2 (111) Figure 7: A plot of surface coverages and partial pressures during a microkinetic analysis of the ZrO 2 surface at 700 K with 1:1 H 2 O:CO ratio and 2 bar s −1 flow rate. For the sake of clarity, intermediates whose coverage was < 0.01 monolayers throughout the entirety of the simulation time were not plotted.
As expected based on experimental results on the bare zirconia, 15 no gas-phase H 2 or CO 2 are produced at any temperature or gas composition studied. At 300 K, the surface is totally covered by molecular water. When the temperature reaches 350 K water dissociation becomes feasible and ∼ 5 % of water molecules dissociate, at 500 K already ∼ 60 % of water dissociates into OH and H. Formate formation begins at 600 K. Figure 7 shows how formate accumulates on zirconia as a function of time at 700 K. Even at the beginning of the simulation, all molecular water has dissociated and the overall coverage of dissociated water fragments is 0.8 ML. The formate coverage gradually increases reaching 0. 6 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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 experimental results, where formate formation was found to occur at the expense of terminal OH, while the concentration of "multicoordinated" OH groups remained unchanged. 15 We note that further increase of temperature allows to reach steady state faster but leads to decrease of formate concentration due to fewer terminal OH groups available.

WGS on Rh(111)
On the planar Rh(111) surface, no formation of either product, CO 2 or H 2 , was observed independent of reaction conditions and gas ratios studied in microkinetic simulations. This is because the surface is fully occupied by CO molecules with a coverage of 1 ML, which prevents water adsorption. Such a high CO coverage was also observed in microkinetic simulations for the WGS on Al 2 O 3 supported Ni nanoparticles. 18 This was attributed to excluding lateral interactions between adsorbates from the model. A high CO coverage was also reported for Pt(111), when lateral interactions were absent. 32 The CO adsorption energy on metal facets is known to be coverage dependent due to the repulsive lateral interactions between the CO molecules. 77 This is also true for rhodium, 78 still, 0.75 ML CO coverage has been reported on rhodium single crystals under ultra high vacuum 79 . To include the coverage dependence of CO adsorption into our model, lateral interactions between CO molecules were taken into account by employing the mean-field approach according to equation 2.0.4 in SI.
Simulations with coverage dependent CO adsorption show that the coverage varies from 0.5 ML to 0.8 ML with temperature ranging from 800 to 300 K, but is independent of the H 2 O:CO ratio. Even though there are now free surface sites available for water adsorption and dissociation, no product formation is observed. Our DFT binding energy for water was calculated at 0.25 ML, but a recent DFT study demonstrates that water adsorption energy on metals becomes more exothermic at higher coverage. 73 We hypothesize that weak water adsorption, rather than poisoning of Rh by CO, prevents the WGS reaction from proceeding.
To unravel the reaction mechanism, we use the degree of rate control (DRC) analysis. 72 The summation of the X RC values of all steps in the microkinetic model should equal to unity, meaning that in the case of a single rate determining step, for that step X RC = 1. As DRC can only be evaluated when products are formed, we analyzed the reaction mechanism for enhanced H 2 O adsorption energy of -0.82 eV, that is 0.5 eV more exothermic than the DFT computed value, but close to the value reported for a high water coverage. 73 This gives a 7 % CO conversion at 700 K, which is substantially lower than at the interface, as shown in the next section. The obtained DRC values are X RC, R2 = 0.418 and X RC, R7 = 0.403 for water dissociation and CO oxidation steps, respectively. These findings show that WGS follows a surface redox mechanism, which is limited by two elementary steps. The negative DRC value of X RC, R3 = −0.241 for the complete water dissociation step indicates that the atomic oxygen needed for CO oxidation must originate from the disproportionation reaction.
Interestingly, according to the DFT results, the carboxyl pathway is energetically less costly on Rh(111) than the surface redox reaction. The DRC value for the carboxyl dissociation step is only X RC, R10 = 0.129, which means it is not as rate controlling as the CO oxidation step. The carboxyl formation step has a small negative DRC value, indicating that this step slightly inhibits rate of H 2 production, which in turn means it cannot be part of the dominating pathway. The reason why the carboxyl reaction is unfavorable is likely due to the low availability of OH groups. After formation, the OH groups are immediately consumed by the disproportionation reaction, which is more favorable than COOH formation, due to the more exothermic reaction energy and zero kinetic barrier.
The CO conversion and coverages of the most abundant surface species, CO, O, and H are plotted against H 2 O adsorption energy in Figure 9. We observe that CO conversion increases with increasing O coverage, while CO coverage remains constant. As water adsorption becomes more exothermic, partial dissociation becomes more feasible. This in turn provides more OH groups for the disproportionation reaction generating the atomic oxygen 25 Figure 9: CO, O, and H coverages, and CO conversion, %X CO , at steady state WGS plotted against H 2 O adsorption energy. The CO adsorption energy was kept at the original DFT value. Note that the lines connecting the data points are for the purpose of guiding the eye only, and do not indicate interpolation of data. The reaction conditions were kept at 700 K with 1:1 H 2 O:CO ratio and 2 bar s −1 flow rate. and leading to higher CO conversion.

WGS at the Rh/ZrO 2 -interface
At the Rh/ZrO 2 -interface, the production of H 2 and CO 2 begins at 600 K, while at lower temperatures CO and OH are present but do not react further. Since CO binds on the Rh top site at the interface, lateral interactions between CO molecules were treated in the same manner as for Rh(111), but we note that the reaction does proceed even without them. Figure   10 illustrates steady state CO conversion,%X CO , plotted against temperature at differing gas compositions. In general, CO conversion increases with an increasing water:gas ratio, and the maximum conversion, 75%, is reached at 700 K with the ratio 10:1. The simulations demonstrate negligible surface coverage for formate between 300 K and 500 K at all water:gas ratios. However, similar to the zirconia support, formate appears at 600 K and its coverage is ca 0.1 ML at all gas ratios, falling below 0.1 ML at 700 K.
Naturally, lowering the dissociation barrier of such a species will have a positive effect on the overall reaction rate. Based on this DRC data, we conclude that at the Rh/ZrO 2 -interface the WGS follows the associative reaction mechanism with carboxyl as the reaction intermediate.  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 high water to gas ratio in Figure 11. At both gas compositions, CO coverage decreases with increasing temperature, while OH coverage goes through a maximum and decreases at high temperature. Higher partial pressure of water leads to a higher coverage of OH species, which in turn facilitates higher CO conversion. At 800 K, the conversion is limited by the low coverage of reactants, in particular OH.
Next, we compare the coverages of the most abundant reactive species on the rhodium and interface models. The reaction conditions were selected to be 700 that CO partial pressure has a minor effect on the overall reaction rate, while the influence of water partial pressure is substantial. Our result for CO agrees very well with a previous experimental result, while the effect of water partial pressure is overestimated in our study. 22 A similarly over estimated H 2 O order was obtained from a WGS microkinetic model on Pt(111). 32 Correctly capturing the negative CO 2 reaction order is typically challenging, 8,40,81 but, in agreement with experiments, 22 we find a negative reaction order for CO 2 , our value −0.12 being less negative than the experimental one. In contrast, we also find a negative order with respect to H 2 , which was previously experimentally determined to be close to zero. 22 Unfortunately the limited amount of experimental kinetic data makes comparison between the predicted and experimental reaction orders difficult, as the values typically 28  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 depend on the reaction conditions, and prevents the use of selective refinement schemes. 8,82,83 Discussion Challenges in modeling catalysis at the metal/oxide-interface In general, construction of computational models for metal/support-interfaces presents a number of challenges for atomistic simulations. One such challenge is the selection of an exchange-correlation functional, which determines the quality of the DFT calculations. It  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 both in the effectively screened "solid regions" and at the unscreened tail of the "molecular regions". Even if the exact DFT functional is available, the evaluation of free-energies is not easy and may contribute to observed differences between experimental and computational results. 96 The construction of an atomistic model for the catalyst/support-interface introduces ad-

Computational insights into WGS on zirconia supported rhodium
The present study provides new insight into the reaction mechanism and active domain of the water-gas reaction over the zirconia supported rhodium catalyst. According to our DFTmicrokinetic analysis, the interface is the most active domain of the Rh/ZrO 2 catalyst, which is not easily elucidated from DFT results alone because the energy differences between the competing pathways are minor. We predict that the support and bulk rhodium are inactive 30  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  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 for WGS. Only when the interface is explicitly included, product formation is observed. The turn over frequency (TOF) for the reaction at the interface is 38 times greater than the TOF over Rh(111) at identical reaction conditions (700 K with 1:1 H 2 O:CO ratio and 2 bar s −1 flow rate), even when using the modified H 2 O adsorption energy for Rh(111). This implies that while the interface sites make up a minority of the total surface sites, they are still responsible for the majority of the WGS activity. We also note that decreasing particle size leads to a greater number of interfacial sites, and in turn results in greater catalytic activity. Conversely, as metal sintering has been previously identified as a possible deactivation mechanism for WGS over Rh/ZrO 2 , 22 this further supports the above interpretation that the interface is the catalytically active domain. Thus, neglecting the oxide-metal interface in computational models might lead to incorrect conclusions on the activity of different catalysts towards WGS. Even if ordering of catalytic activity of different metals can be correctly captured with simplified models, neglecting the metal-support interface even on nominally inert supports may lead to unsatisfactory agreement between computational and experimental results. As observed in this work and previously, 8,9,18,19,51 the interfacial sites can be paramount for the overall catalytic behavior. As such, over-simplified computational models might need refinement to predict correct chemistry and to identify promising bifunctional materials not captured by existing scaling relations. 20,52 As shown in the present work, and previously experimentally, 15  to increase the CO conversion by less than 1 % for the conditions and gas compositions considered. 22 We note that higher content of WGS products in the inlet gas substantially enhances the methane yield. 7 The slightly lower CO conversion compared to the measured value could be due to approximations and models employed in calculations, but also due to 32 the differences in reaction conditions, catalyst structure, and reactor model.

Conclusions
The first-principles microkinetic model presented herein shows how the water-gas shift (WGS) reaction proceeds over a ZrO 2 -supported Rh catalyst. We find that the metal-support interface is the most active domain of the Rh/ZrO 2 system giving highest CO conversion to hydrogen and carbon dioxide. The reaction proceeds along the associative reaction mechanism via a carboxyl intermediate. Interestingly, the carboxyl coverage remains negligible under all reaction conditions and gas ratios studied since it reacts onward immediately, explaining why carboxyl isn't readily observed in experiments.
DFT-computed CO adsorption energy on Rh(111), representing the facets of Rh nanoparticles, leads to 1 ML CO coverage in the microkinetic analysis and renders the surface unable to run the reaction. Including mean-field lateral interactions between CO molecules reduces CO coverage but does not initiate WGS reaction, which is only achieved if the exothermicity of water adsorption is enhanced. In contrast to the interface, the reaction proceeds now via the surface redox route with water activation and CO oxidation as the rate controlling steps.
However, even with this modified water adsorption energy, the CO conversion on Rh (111) remains much lower than that achieved at the interface.
The exposure of m-ZrO 2 (111) to reactant gases leads to the formation and accumulation of formate from 600 K onwards in microkinetic simulations. Formate on the bulk zirconia is shown to act as a spectator species. Calculations show that formate diffusion to the metalsupport interface and its subsequent decomposition are feasible. Even so, the rate of reaction for this process is expected to be slow, and therefore it provides only a minor contribution towards the WGS products, as compared to the carboxyl pathway at the interface.
Our computational study provides an atomically refined model of WGS on zirconia supported rhodium nanoparticles. We show that different parts of the system have distinctly 33 different roles in the reaction. Although the bulk oxide support and metal particles do not directly take part in the WGS reaction, it is fair to say that they both have a major effect on the reaction mechanism through the formation of the metal-oxide interface. Based on the obtained atomistic information from the combined DFT and microkinetic model, surface coverages and CO conversions were calculated and related to the elementary steps. The obtained results explain and resolve experimental observations and interpretations related to the formation and role of formate as well as the absence of carboxyl. This was achieved by including the support, the metal catalyst, and the interface within the computational model.
The present DFT and microkinetic study further strengthens the argument that design of bifunctional catalysts must not rely on simplified models for activity and selectivity as this can potentially lead to the exclusion of prominent material combinations and inaccurate designing principles for a multitude of catalytic reactions.