Hydrodeoxygenation of guaiacol over orthorhombic molybdenum carbide: a DFT and microkinetic study

Abstract

The hydrodeoxygenation of guaiacol is modelled over a (100) β-Mo₂C surface using density functional theory and microkinetic simulations. The thermochemistry of the process shows that the demethoxylation of the guaiacol, to form phenol, will be the initial steps, with a reaction energy of 29 kJ mol⁻¹ (i.e. endothermic) and a highest activation barrier of 112 kJ mol⁻¹. Subsequently, the dehydroxylation of the phenol, which has a rate-determining activation barrier of 145 kJ mol⁻¹, will lead to the formation of benzene, with an overall reaction energy for conversion from guaiacol of −91 kJ mol⁻¹ (i.e. exothermic). The sp² and sp hybridized carbon atoms of the molecular functional groups are found to dissociate on the surface with minimum energy barriers, while the hydrogenation of the adsorbed molecules requires higher energy. The microkinetic modelling, which is performed considering typical reaction conditions of 500 to 700 K, and a partial pressure ratio for H₂ : guaiacol of 1, shows quick formation and accumulation of phenol on the surface with increasing temperature, although high temperatures mitigate the guaiacol adsorption step. Based on simulated temperature programmed desorption (TPD), maximum conversion of guaiacol can be expected at 70% surface coverage of this species.

---

1. Introduction

In the quest for a viable renewable source of energy, biomass feedstocks have emerged as a strong contender because of the abundance of natural flora on the planet. The bio-oil obtained from the biomass, when properly treated, can be used as a direct substitute to fossil fuel;^1,2 however, the cost of such treatment is currently preventive for commercial application. In particular, bio-oil obtained from the thermal treatment of biomass contains oxygen carrying compounds like guaiacol, anisole, and ferulic acid, and these oxy-compounds need to be reduced before being suitable as fuel.³ The reduction step is necessary because the oxygen carrying compounds are known to provide unfavorable properties to the crude bio-oil (such as low pH, high viscosity, etc.), and reduce its energy density.⁴

The reduction of the oxy-compounds can be achieved through a process called hydrodeoxygenation (HDO), where molecular hydrogen is reacted at moderate temperature (300–500 °C) and high pressure (up to 30 bar) with the bio-oil,⁵ and oxygen is subsequently removed in the form of water. Unfortunately, the conditions outlined for HDO are economically prohibitive for scale-up;⁶ thus, catalysts are being actively sought to facilitate the reaction at milder conditions. Several different classes of metals have been tested as catalysts for the HDO process. While most of the transition metals are poor catalysts, due to low yield and high rates of char formation, the noble metals Ru, Pt and Pd have shown good performance;⁷ however, their low earth abundance and subsequent high cost prevents these materials being considered as commercial catalyst towards HDO.

To overcome the challenge of identifying effective and affordable catalysts, scientists are now considering cheaper materials with catalytically relevant characteristics similar to those of the noble metals. Nørskov et al.⁸ have shown that the performance of a catalyst is significantly influenced by the electronic structure at its surface, with the position of the d-band center, relative to the Fermi level, identified as a suitable indicator of the materials adsorption strength and catalytic activity.⁹ Such direct relation between the catalytic performance and the d-band position of the transition metal is due to the bonding and anti-bonding states formed when the adsorbate approaches the surface, where hybridization between the frontier orbitals of the adsorbate and s and d-orbital of the metal occurs; the adsorption strength is therefore dependent on the relative position of the metal states with respect to the Fermi level. The d-band theory has been quite successful in explaining, for instance, the trends observed in the oxygen reduction reactions (ORR),^10–13 and the d-band center of catalytic materials has also been used as a descriptor for property prediction in transition-metal catalysed HDO studies.^14–17 The transferability of the d-band approach is now of consideration for alternative composite catalysts, such as carbides and nitrides. As an example, molybdenum carbide (Mo₂C) has similar d-band characteristics to the catalytically active noble metals (Ru, Pt and Pd), due to the carburization of molybdenum;¹⁸ specifically, the introduction of carbon between the Mo atoms increases their separation, which shifts the d-band center closer to the Fermi level, akin to the properties of the highlighted noble metals.

Considering the valorization of biomass, the crude bio-oil obtained after fast pyrolysis has a relatively high concentration of guaiacol (5–15% of the phenolic fraction).² Guaiacol also contains both a methoxy and a hydroxyl functional group, which represents most of the oxygen-containing functional groups in bio-oils. In combination, these factors make guaiacol a convenient model molecule for studying the HDO of the bio-oil phenolic fraction. The upgrading of guaiacol over various catalysts, including the noble metals, and materials containing molybdenum and carbon with varying stoichiometries, has been considered in experimentation by several groups.^19–28 Ma et al.²⁹ upgraded guaiacol over α-MoC supported on activated carbon, with combined selectivity to phenols and alkylphenols of over 85%, and conversion at 87%. Similarly, Jongerius et al. obtained a combined selectivity of up to 87%, and a conversion greater than 99%, for the HDO of guaiacol over carbon nanofiber supported Mo₂C;³⁰ in this latter case, the upgrading was conducted at mild conditions of 5 bar of hydrogen pressure and temperatures of 300 to 375 °C. The Mo₂C catalyst outperformed a tungsten carbide (WC) system with respect to yield of completely deoxygenated products, such as benzene. Indeed, Siaj et al.³¹ found acetaldehyde to dissociate selectively from the carbonyl bond over β-Mo₂C surface in their experiments, with the mechanism explained subsequently by Martínez et al.³² Moreira et al.¹⁸ studied the efficacy of a Mo₂C catalysts, supported again on carbon nanofiber, during the HDO of guaiacol, with phenol and cresol observed as the main products; the formation of phenol is proposed to have been initiated by the demethylation of guaiacol.

Experiments have clearly shown that Mo₂C is a potentially effective and affordable catalyst for HDO reactions; however, the exact mechanism and the rate controlling parameters of the HDO process over the Mo₂C surfaces remain unknown despite being necessary for catalyst design. Therefore, in this study, the reaction pathway for upgrading guaiacol over the Mo₂C surface is investigated using density functional theory. A network of routes is proposed for guaiacol conversion to benzene, with important intermediates such as phenol and catechol included. The thermochemistry of all the proposed reactions is presented (Section 3) and subsequently, microkinetic modelling is applied to determine the kinetic parameters (Section 4) of the upgrading process.

2. Computational details

The orthorhombic (β) structure of molybdenum carbide was selected as base structure for our studies, as the excellent thermal stability of the material³³ makes it suitable for catalytic applications over a wide temperature range. The (100) surface with a molybdenum termination shows favourable hydrogen adsorption properties,³⁴ and has been reported as a suitable surface for HDO reactions in experimental works,^29,35 and thus is considered herein. The Mo terminated (100) surface also has the highest surface metal density (0.130 atom per Ų) of all terminations of the (100) surface, resulting in a high coordination of surface Mo atoms. As a result, the electron-state fluctuations, which are involved in the breaking and formation of chemical bonds during catalytic reactions and are responsible for high turnover, are the highest at the Mo terminated surface.³⁶ To simulate the β-Mo₂C (100) surface, a 4-layer slab model with a 3 × 3 supercell was considered in a periodic environment (Fig. 1).

Fig. 1 (a) Top view of the β-Mo₂C (100) slab showing adsorption sites, (b) side view of the β-Mo₂C (100) slab showing the vacuum, and (c) methoxy group, hydroxyl group and the aromatic ring of guaiacol; where the grey, teal, red and white coloured atoms represent carbon, molybdenum, oxygen and hydrogen, respectively.

All energy calculations were performed under the density functional theory (DFT) framework using the “Fritz Haber Institute ab initio molecular simulations” (FHI-aims) software package³⁷ in combination with the “Atomic Simulation Environment” (ASE) Python package³⁸ for geometry handling. After convergence testing with respect to surface energies, constraints were deemed appropriate for the bottom two layers of the slab, thus maintaining the long-range bulk structure, whilst the top two layers were unconstrained during all adsorption and reaction modelling. A 10 Å vacuum was added above and below the plane of the slab (i.e. total vacuum of 20 Å), which prevents spurious self-interaction errors. To counterbalance any dipole arising out of the imbalanced electric field at the surface, a dipole correction was also applied. A converged 5 × 5 × 1 k-grid was applied for the periodic condition calculations. The PBE³⁹ functional with the Tkatchenko-Scheffler van der Waals correction⁴⁰ was used for electronic structure calculations, along with the “light” basis set (version: 2010) of the FHI-aims package.⁴¹ In addition, zeroth order regular approximation (ZORA)³⁷ scalar corrections were incorporated to account for relativistic effects. The choice of all these parameters is a result of a systematic study described previously.⁴²

Geometry optimization of all the structures was conducted using the trust region method⁴³ until the force on each atom was less than 0.01 eV Å⁻¹. For the transition state calculations, a minimum of 7 images were used in the nudged elastic band (NEB) method calculations⁴⁴ with the molecular dynamics based fast inertial relaxation engine (FIRE) optimization algorithm.⁴⁵ The convergence criteria for the NEB calculations was set to a requirement of forces on each atom being less than 0.05 eV Å⁻¹. This was followed by finite-difference frequency calculations with 0.01 Å displacement of atoms for the transition state structures to confirm their validity, with one imaginary frequency confirming a first-order saddle point. For cases where none or multiple imaginary frequencies were obtained, a more exhaustive machine-learning NEB method⁴⁶ was used with 15 images to determine the correct transition state structure. In specific cases, where obtaining the transition state was challenging, a complementary approach was used: the distance, d, between atom A (or a group of atoms) that reacts with another atom B (or a group of atoms) was calculated, and then divided into n equal parts (n =20 in this case). Atom(/group) A is then placed at a distance of d/n from atom(/group) B, for all values of n, generating a pathway that contains n different geometry samples, and then structural optimization is performed for each geometry, with a cut-off force of 0.05 eV Å⁻¹ on each atom, with the distance between A and B constrained. On plotting energy versus d for all the optimized geometries, the energy profile can then be used to identify the transition state structure. Once the transition state structure is obtained, again finite-difference frequency calculations were applied to confirm the nature of the first-order saddle point through identification of a single imaginary frequency.

The adsorption energy (E_ads) of appropriate reaction processes was calculated as:

E_ads = E_{Mo2C+Molecule} − E_Mo2C − E_Molecule (1)

where E_Molecule is the energy of the adsorbate molecule in the gas phase, E_Mo2C is the energy of the bare slab surface, and E_{Mo2C+Molecule} is the energy of the combined system where the adsorbate is on the slab surface. For the reactions describing the desorption of molecules from the surface, the desorption energy (E_des) was calculated as E_des = −E_ads.

The microkinetic modelling was conducted using our in-house code,⁴⁷ with details outlined elsewhere.⁴⁸ To describe briefly the details here: translations, rotations and vibrations were used to calculate the thermodynamic parameters, such as entropy, enthalpy and Gibbs free energy, based on a statistical thermodynamics approach.⁴⁹ The derived energetic quantities were then used to calculate the rate constants for all the reactions using the transition-state theory (TST) approximation of Eyring, Evans and Polanyi^50,51 as shown in eqn (2).

Here, A₀ is the pre-exponential factor, ΔG^‡ is the activation free energy of the reaction, k_B is the Boltzmann constant, h is Planck's constant, and Q_TS and Q_r are the partition functions of the transition state and reactant, respectively. The microkinetic model builds on transition state theory with an improved description of tunnelling barriers.⁴⁸ It assumes that each site on the surface is identical, that the adsorbed species are adsorbed randomly on the surface and do not interact laterally, and that each reaction considered in the model is an elementary reaction unhindered by any mass transfer and heat transfer resistance.⁴⁹ For adsorption reactions, the rate constants were calculated using the Hertz–Knudsen relation.⁵² Finally, the rate of reaction was described for each individual step and the system of ordinary differential equations (ODEs) was solved to obtain a steady state solution.

3. Results and discussions

The mechanism proposed herein initiates with the adsorption of guaiacol (GUA) on the catalyst surface. Guaiacol was considered interacting with the surface in three different orientations via: the aromatic ring, the methoxy group, or the hydroxyl group, as shown in Fig. 1(c). The adsorption was conducted on 5 different sites on the surface (atop, bridge, fcc, Mo-hcp and C-hcp, as illustrated in Fig. 1(a)), with the adsorbed structures and energetics shown in Tables S1–S3 of the ESI.†

The strongest adsorption of guaiacol is via the aromatic ring over the C-hcp position of the surface, with an adsorption energy of −4.67 eV, which is at least 2.21 eV stronger than reported for the precious metal surfaces such as Ru (E_ads = −2.46 eV),⁵³ Pd (E_ads = −1.43 eV)⁵⁴ and Pt (E_ads = −2.41 eV);⁵⁵ the adsorption energy for guaiacol on tungsten carbide is closer to our calculated results, having been previously reported as −3.04 eV,⁵⁶ though we note this is still 1.63 eV weaker than the molybdenum carbide. The difference in the adsorption energy could be a result of use of different functionals in different studies. In several cases during our investigation of orientations for guaiacol adsorption, such as interaction in the atop position with the methoxy group, and in the bridge position via the aromatic ring, the adsorbed molecule rearranges in the optimization process to a position over the C-hcp position and interacting via the aromatic ring. The preference of the hollow site (similar to the C-hcp site) for adsorption of aromatic species, including guaiacol, is well documented.⁵⁷ As the C-hcp position was confirmed as the most stable adsorption site, all the subsequent calculations considering conversion of guaiacol were conducted at this position.

In order to validate the observations, calculations were conducted to obtain the sticking probability (σ(T)) of aromatic compounds as:

where q^reac_2D-vib is the 2D vibrational partition function of the adsorbate and q^cell_2D-tran is the 2D translational partition function of the surface, and are obtained from the frequencies of the system.⁴⁹ For guaiacol, our kinetic rate calculations give σ as 0.76 at 170 K, which decreases to the orders of 10⁻⁴ and 10⁻⁷ at 300 K and 500 K, respectively. Thus, the rate constant for adsorption decreases with increase in temperature; or, alternatively, low temperatures are favorable to promote the adsorption process. Whilst comparable literature for guaiacol is unavailable, σ for phenol is very close to unity at 90 K over Pt (111) and at 150 K over Ni (111),⁵⁸ which agrees with the temperature of 160 K below which unity is observed in our calculations for guaiacol. In contrast, over a Ag (111) surface, σ for phenol is 0.56 at 163 K,⁵⁹ which is 0.20 lower than our sticking coefficient at a similar temperature (0.76 at 170 K).

3.1. Energy profile of the upgrading routes

The HDO of guaiacol has been studied extensively over monometallic^19,22,24,60 and bimetallic surfaces.^21,27,61,62 Based on the reported literature, the possible pathway for upgrading could be guaiacol → catechol → phenol → benzene. Therefore, the proposed mechanism for hydrodeoxygenation involves the hydrogenation and subsequent deoxygenation of guaiacol in a series of elementary steps. Full saturation of the aromatic ring with hydrogen is reported as feasible at high hydrogen pressure and low temperature (<573 K) over noble metals;²⁸ however, since completely saturated products are not desirable, and observed outside the standard HDO conditions, we focus on the upgrading of the adsorbed guaiacol beginning in two different ways: hydrogenation at the α- or the β-position of the aromatic ring (Fig. 2, structures 2 and 3), or dissociation of the the methoxy and hydroxyl groups from the molecule. In combination, there are therefore seven routes (R1–R7) by which the upgrading can initiate (Fig. 2).

Fig. 2 Reaction scheme of all the elementary steps for guaiacol upgrading; numbers given in red are the forward kinetic barriers of each reaction, in eV, whilst blue and black text show the structure and reaction numbers, respectively. Equivalent structures are highlighted with identically coloured backgrounds.

With respect to kinetic barriers, the dissociation of H from the hydroxyl group (R7) is the least energy demanding (0.62 eV), whereas thermodynamically, R6 is the most favourable reaction (ΔE = −2.29 eV). The deprotonation of the hydroxyl group was previously observed as having the lowest energy barrier among R1–R7 on noble metal surfaces,^28,53 showing similarity to our results; the reported barrier over Pt (111) surface is 0.37 eV,²⁸ and it is 0.29 eV over the Ru (0001) surface⁵³ for the same reaction. The next most accessible reaction in our calculations is the cleavage of the methoxy group (R4), which has an activation energy (E_a) of 0.90 eV. The kinetic rate calculations show the rate constant for R7 to be 100 times faster than R4, suggesting faster production of structure 8 (6-methoxycyclohexa-2,4-dienone) over structure 5. The direct hydrogenation of the molecule via R1 and R2 is highly energy demanding with a kinetic barrier of 1.15 eV and 2.59 eV, respectively. Since the hydrogen cannot be transferred to the α- or β-position atom directly from the surface, due to steric hindrance, a concerted reaction mechanism is considered in these reactions. In the concerted process, a surface hydrogen migrates to the aromatic ring, and simultaneously a hydrogen from the ring shifts to the α- or β-position. Initially, the high activation barrier was suspected as arising due to hydrogen diffusion on the surface; however, R49 (Table 1), which describes the diffusion and association of hydrogen atoms to form molecular hydrogen on the surface, shows that the hydrogen diffusion is barrierless on the β-Mo₂C(100) surface. Therefore, the activation energy observed is associated with the second step in the concerted process, i.e., the migration of H from the aromatic ring to the α- or β-position.

Table 1 Activation energy (E_a), reaction energy (ΔE), pre-exponential factor (A₀), forward rate constant (k_f) and activation free energy (ΔG^‡) for reactions described in the reaction scheme. The notation used for the reactants (Xr), transition states (TSX) and products (Xp) uses X to represent the reaction number. The surface is denoted as * in the reaction process, such that e.g. Xr* denotes an adsorbed reactant. The structures, as shown in Fig. 2, are denoted as SY, where Y represents the structure number

Rxn. No	Reactions	E a (eV)	ΔE (eV)	500 K			600 K			700 K
				A 0	k f (s^−1)	ΔG^‡ (eV)	A 0	k f (s^−1)	ΔG^‡ (eV)	A 0	k f (s^−1)	ΔG^‡ (eV)
a Reactions with transition state calculated by the complementary approach to NEB.
Ads	GUA + * ≫ GUA*	—	−4.67	1.30 × 10^4	4.19 × 10^−3	−0.05	1.18 × 10^4	4.46 × 10^−4	−0.04	1.10 × 10^4	6.45 × 10^−5	−0.04
R1	1r* > TS1 > 1p*	1.15	0.64	4.47 × 10^13	1.78 × 10^5	0.26	4.46 × 10^13	2.15 × 10^5	0.25	4.39 × 10^13	2.50 × 10^5	0.25
R2	2r > TS2 > 2p*	2.59	0.86	8.18 × 10^13	1.78 × 10^−4	0.33	9.37 × 10^13	2.35 × 10^−4	0.33	1.05 × 10^14	3.13 × 10^−4	0.33
R3^a	GUA* > TS3 > 3p*	1.04	−0.69	4.03 × 10^13	4.82 × 10^6	0.96	4.69 × 10^13	6.33 × 10^6	0.95	5.34 × 10^13	8.36 × 10^6	0.95
R4	GUA* > TS4 > 4p*	0.90	0.59	7.62 × 10^13	9.65 × 10^7	0.82	8.67 × 10^13	1.32 × 10^8	0.81	9.60 × 10^13	1.79 × 10^8	0.80
R5	GUA* > TS5 > 5p*	1.02	−0.49	2.26 × 10^13	6.76 × 10^6	0.91	2.47 × 10^13	7.09 × 10^6	0.91	2.70 × 10^13	7.78 × 10^6	0.91
R6	GUA* > TS6 > 6p*	1.44	−2.29	1.41 × 10^15	1.33 × 10^6	1.25	2.09 × 10^15	3.51 × 10^6	1.22	2.87 × 10^15	9.13 × 10^6	1.18
R7	GUA* > TS7 > 7p*	0.62	−0.65	1.21 × 10^13	2.01 × 10^9	0.52	1.13 × 10^13	1.56 × 10^9	0.54	1.07 × 10^13	1.24 × 10^9	0.55
R8	1p* > TS8 > 8p*	0.75	−0.46	2.83 × 10^13	4.32 × 10^8	0.67	2.93 × 10^13	4.59 × 10^8	0.67	3.03 × 10^13	4.98 × 10^8	0.66
R9	1p* > TS9 > 9p*	0.37	−2.05	2.77 × 10^14	3.68 × 10^12	0.26	3.31 × 10^14	6.24 × 10^12	0.24	3.79 × 10^14	1.04 × 10^13	0.22
R10	2p* ≫ 10p*	—	−2.16	2.68 × 10^13	9.58 × 10^28	−2.16	3.04 × 10^13	1.23 × 10^29	−2.17	3.37 × 10^13	1.58 × 10^29	−2.18
R11	11r* > TS11 >11p*	0.97	−0.29	7.06 × 10^13	2.13 × 10^7	−0.03	8.24 × 10^13	3.07 × 10^7	−0.03	9.31 × 10^13	4.33 × 10^7	−0.04
R12	12r* > TS12 > S26*	1.16	−0.29	6.48 × 10^13	7.57 × 10^5	−0.02	7.50 × 10^13	1.08 × 10^6	−0.02	8.38 × 10^13	1.49 × 10^6	−0.03
R13	S6* > TS13 > 13p*	0.79	−0.22	1.86 × 10^15	8.08 × 10^8	0.88	4.53 × 10^15	9.85 × 10^9	0.79	1.04 × 10^16	1.19 × 10^11	0.69
R14	S6* > TS14 > 14p*	0.82	−1.58	5.70 × 10^14	3.63 × 10^8	0.86	1.24 × 10^15	3.05 × 10^9	0.78	2.64 × 10^15	2.63 × 10^10	0.69
R15	S6* > TS15 > 15p*	0.49	−2.08	5.63 × 10^14	2.43 × 10^10	0.61	1.20 × 10^15	2.16 × 10^11	0.52	2.48 × 10^15	1.92 × 10^12	0.43
R16	16r* > TS16 > 16p*	1.81	1.43	1.94 × 10^14	8.80 × 10^1	−0.03	2.30 × 10^14	1.43 × 10^2	−0.03	2.61 × 10^14	2.27 × 10^2	−0.03
R17	17r* > TS17 > 17p*	2.46	0.52	1.32 × 10^13	1.11 × 10^−4	0.58	1.33 × 10^13	9.74 × 10^−5	0.57	1.33 × 10^13	8.60 × 10^−5	0.57
R18	S7* > TS18 > 18p*	1.03	−1.19	8.30 × 10^13	1.44 × 10^7	0.94	9.39 × 10^13	2.00 × 10^7	0.93	1.04 × 10^14	2.80 × 10^7	0.91
R19	S8* > TS19 > 19p*	0.39	−1.99	7.42 × 10^13	5.60 × 10^11	0.29	8.75 × 10^13	7.82 × 10^11	0.28	1.00 × 10^14	1.09 × 10^12	0.27
R20	S8* > TS20 > 20p*	0.78	−0.54	3.66 × 10^13	7.88 × 10^8	0.65	4.12 × 10^13	9.15 × 10^8	0.65	4.59 × 10^13	1.10 × 10^9	0.64
R21	S9* ≫ 21p*	—	−1.94	1.06 × 10^13	2.45 × 10^27	−1.99	1.15 × 10^13	2.49 × 10^27	−1.99	1.23 × 10^13	2.57 × 10^27	−1.99
R22	S11* > TS22 > 22p*	1.02	−0.52	2.57 × 10^13	8.67 × 10^6	0.90	2.85 × 10^13	9.34 × 10^6	0.90	3.14 × 10^13	1.05 × 10^7	0.90
R23	S18* ≫ 23p*	—	−2.05	2.00 × 10^13	2.97 × 10^28	−2.11	2.18 × 10^13	3.33 × 10^28	−2.11	2.35 × 10^13	3.73 × 10^28	−2.11
R24	16p* > TS24 > 24p*	0.18	−1.12	6.01 × 10^13	9.12 × 10^12	0.11	6.83 × 10^13	1.22 × 10^13	0.10	7.60 × 10^13	1.65 × 10^13	0.09
R25	25r* ≫ 25p*	—	2.98	1.58 × 10^13	1.73 × 10^−8	0.47	1.55 × 10^13	1.54 × 10^−8	0.47	1.52 × 10^13	1.37 × 10^−8	0.46
R26	17p* > TS26 > 26p*	0.09	−1.87	1.64 × 10^13	4.40 × 10^12	0.08	1.68 × 10^13	4.56 × 10^12	0.08	1.71 × 10^13	4.73 × 10^12	0.08
R27	27r* > TS27 > S19*	1.15	0.23	7.30 × 10^13	9.72 × 10^5	−0.65	8.63 × 10^13	1.44 × 10^6	−0.65	9.81 × 10^13	2.07 × 10^6	−0.66
R28	S12* ≫ 28p*	—	−2.31	1.87 × 10^13	2.54 × 10^30	−2.38	2.09 × 10^13	2.86 × 10^30	−2.38	2.30 × 10^13	3.23 × 10^30	−2.38
R29	29r* > TS29 > 25p*	1.22	3.61	1.07 × 10^13	5.05 × 10^5	0.24	1.20 × 10^14	7.15 × 10^5	0.24	1.30 × 10^14	9.89 × 10^5	0.23
R30	30r* > TS30 > S26*	2.19	1.41	5.97 × 10^13	3.00 × 10^−2	−0.02	6.86 × 10^13	4.11 × 10^−2	−0.03	7.64 × 10^13	5.54 × 10^−2	−0.03
R31	25p* ≫ 31p*	—	−4.47	1.96 × 10^13	4.79 × 10^44	−4.36	2.19 × 10^13	6.78 × 10^44	−4.37	2.40 × 10^13	9.44 × 10^44	−4.39
R32	32r* > TS32 > S7*	2.28	1.69	1.89 × 10^14	2.99 × 10^−2	0.00	2.25 × 10^14	4.95 × 10^−2	−0.01	2.58 × 10^14	7.98 × 10^−2	−0.01
R33^a	S26* > TS33 > 33p*	1.50	−0.73	7.27 × 10^13	6.05 × 10^3	1.40	8.62 × 10^13	8.53 × 10^3	1.39	9.91 × 10^13	1.20 × 10^4	1.38
R34	34r* > TS34 > 34p*	2.06	0.57	1.66 × 10^14	5.05 × 10	0.89	1.98 × 10^14	7.46 × 10	0.88	2.29 × 10^14	1.12 × 10^1	0.87
R35	35r* > TS35 > 35p*	0.75	−0.5	5.26 × 10^13	5.33 × 10^8	0.46	6.09 × 10^13	7.34 × 10^8	0.45	6.81 × 10^13	9.88 × 10^8	0.45
R36	34p* > TS36 >36p*	0.5	−1.94	1.83 × 10^13	6.63 × 10^9	0.48	1.84 × 10^13	6.39 × 10^9	0.48	1.84 × 10^13	6.17 × 10^9	0.48
R37	37r* > TS37 >	0.72	0.14	1.90 × 10^13	1.08 × 10^8	0.27	1.92 × 10^13	1.16 × 10^8	0.27	1.93 × 10^13	1.22 × 10^8	0.26
R38	38r* > TS38 >	1.56	1.25	7.62 × 10^13	7.04 × 10^2	0.50	8.21 × 10^13	9.37 × 10^2	0.49	8.62 × 10^13	1.21 × 10^3	0.48
R39	39r* > TS39 > H2O*	2.15	1.65	9.33 × 10^13	1.10 × 10^−1	0.02	1.08 × 10^14	1.58 × 10^−1	0.02	1.21 × 10^14	2.25 × 10^−1	0.01
R40	40r* > TS40 >	0.85	0.17	5.36 × 10^13	3.87 × 10^7	0.12	5.91 × 10^13	5.33 × 10^7	0.11	6.32 × 10^13	7.04 × 10^7	0.10
R41	41r* > TS41 > CH3OH*	2.45	2.14	3.77 × 10^14	6.27 × 10^−3	0.04	4.75 × 10^14	1.17 × 10^−2	0.00	5.66 × 10^14	2.14 × 10^−2	−0.04
R42	S26* ≫ phenol + *	—	4.45	8.91 × 10^24	3.49 × 10^1	3.25	8.24 × 10^24	2.25 × 10^3	3.00	7.41 × 10^24	1.43 × 10^5	2.74
R43	16p* ≫ catechol + *	—	4.26	9.24 × 10^24	6.63 × 10^2	3.08	8.24 × 10^24	4.40 × 10^4	2.82	7.20 × 10^24	2.87 × 10^6	2.56
R44	S29* ≫ benzene + *	—	4.06	2.18 × 10^22	1.18 × 10	3.09	1.79 × 10^22	2.54 × 10^1	2.90	1.46 × 10^22	5.49 × 10^2	2.70
R45	≫ CH4 + *	—	0.38	1.19 × 10^19	5.16 × 10^21	−0.37	1.04 × 10^19	4.18 × 10^22	−0.50	9.06 × 10^18	3.51 × 10^23	−0.64
R46	CH3OH* ≫ CH3OH + *	—	1.18	1.95 × 10^22	1.31 × 10^21	0.16	1.98 × 10^22	3.82 × 10^22	−0.04	1.94 × 10^22	1.13 × 10^24	−0.25
R47	H2O* ≫ H2O + *	—	1.07	4.19 × 10^20	6.31 × 10^18	0.25	4.73 × 10^20	1.22 × 10^20	0.08	5.08 × 10^20	2.33 × 10^21	−0.09
R48	≫ H2 + *	—	0.79	4.91 × 10^18	5.82 × 10^16	0.27	6.88 × 10^18	8.09 × 10^17	0.13	8.82 × 10^18	1.08 × 10^19	−0.01
R49	H* + H* > TS49 >	—	0.00	1.57 × 10^13	1.38 × 10^13	0.01	1.64 × 10^13	1.55 × 10^13	0.00	1.72 × 10^13	1.76 × 10^13	0.00
R50	s19* > TS50 > 37p*	1.59	−0.89	8.27 × 10^13	1.21 × 10^3	1.51	9.67 × 10^13	1.71 × 10^3	1.49	1.10 × 10^14	2.41 × 10^3	1.48
R51	61r* > TS51 > OH*	1.97	1.16	3.77 × 10^13	4.39 × 10^−1	1.94	4.23 × 10^13	5.61 × 10^−1	1.93	4.61 × 10^13	7.02 × 10^−1	1.92

---

The methoxy group of structure 8 undergoes further dehydrogenation in R20, with a reaction barrier of 0.78 eV, to give structure 16; however, structure 16 decomposes rapidly to structure 23 (1,2-benzoquinone) via R28 in a barrierless exothermic step with reaction energy of −2.31 eV. Kinetic rate simulations give a rate constant of ∼10³⁰ s⁻¹, suggesting that the decomposition reaction is very fast. The carbene (CH₂) that separates from structure 16 in R28 sits over a fcc site on the surface with the lowest carbon coordination (i.e., zero); Nagai et al.⁶³ reported similar behaviour for carbon monoxide (CO), with the strongest CO adsorption occurring at the carbon deficient site (i.e., no neighbouring C atoms in the subsurface layer) on the β-Mo₂C surface; experimental studies also report the surface to be selective to the cleavage of the C–O bond.^31,64,65

An oxygen atom of structure 23 (1,2-benzoquinone) can be hydrogenated to produce structure 28 via R32, though a high energy barrier is calculated (E_a = 2.28 eV) and the rate constant for this conversion is low (∼10⁻² s⁻¹). Structure 28 can be formed from two more pathways: (i) the methyl radical from the guaiacol can directly cleave (R6) to form structure 7 (ketophenol) (which is equivalent to structure 28); or (ii) the methoxy group of guaiacol can lose a hydrogen (R5), to give structure 6, followed by the dissociation of a CH₂ moiety (R15) to again produce structure 7/28 (ketophenol). The conversion via R6 is 0.42 eV more energy demanding than R5, and kinetically slower than the R5 and R15 (Table 1). Therefore, the formation of structure 7/28 (ketophenol) will occur predominantly via R5 and R15. The preference of dehydrogenation instead of deoxygenation in the methoxy group of guaiacol has also been reported by Lee et al.⁵⁵ over a Pt (111) surface, with E_a = 0.75 eV for R5 with the help of Brønsted–Evans–Polanyi (BEP) correlation.

From structure 7/28 (ketophenol), there can be three possible reduction routes: in R16, hydrogenation yields structure 13, catechol, with an activation energy of 1.81 eV; in R17, the aromatic ring is hydrogenated with an activation energy of 2.46 eV; and in R18, the hydroxyl group is cleaved in a reaction with activation barrier of 1.03 eV. The kinetic rate modelling returns rate constants for R16 of ∼10¹ s⁻¹ and for R18 of ∼10⁷ s⁻¹, which suggests that the formation of catechol from guaiacol is very slow; furthermore, the high barrier for R17 makes it unlikely. Thus, the reaction process can be concluded as proceeding via R27, with an energy barrier of 1.15 eV, to form structure 19/22, a phenoxyl radical.

The phenoxyl radical (structure 19/22) can also be formed via R2 and R3; the hydrogenation of the β-carbon of guaiacol (R2) has an activation energy of 2.59 eV, and is kinetically very slow (∼10⁻⁴ s⁻¹), whereas R3 (condensation) has a barrier of only 1.04 eV and is kinetically faster. Thus, structure 4 will be formed and converted to structure 11, anisole, via re-hydrogenation of the aromatic ring (R11). Further dehydrogenation of the methoxy group (R22) requires only 1.02 eV of activation energy, which results in an unstable structure 18 that degrades to phenoxyl in a barrier-less exothermic step (ΔE = −2.05 eV).

From the cyclohexadienone, hydrogenation (R30) can yield phenol (Structure 10/26) in a very slow and energy demanding step (E_a = 2.19 eV). The prohibitively high reaction barrier in R30 is in contradiction to experimental studies, which report high selectivity and quick appearance of phenol over pure and supported molybdenum carbides.^18,60,66 Therefore, phenol most likely forms by alternative route(s), of which two routes have been considered herein. One possibility is the direct methoxylation of guaiacol (R4, E_a = 0.90 eV) to give structure 5/20, which is then followed by hydrogenation (R12, E_a = 1.16 eV). The other possibility considered is the hydrogenation of guaiacol (R1) to form structure 2, which has an activation energy of 1.15 eV. The methoxy group can then directly cleave from structure 2 via R9, with a low activation energy of 0.37 eV, producing phenol. R9 is kinetically very favorable, with a rate constant of 3.68 × 10¹² s⁻¹ calculated for 500 K; however, structure 2 can also convert to structure 9 through an activation energy of 0.75 eV, with the unstable CH₂ moiety in structure 9 leading to exothermic degradation (R21) to structure 17. The formation of structure 17 via R21 is kinetically favorable (2.45 × 10²⁷ s⁻¹ at 500 K) as the CH₂ radical dissociates, with a thermodynamic energy change of −1.94 eV. Hydrogenation of structure 17 to structure 21/24 (R29, E_a = 1.22 eV) results in a molecule that decomposes in a highly exothermic step (R31, ΔE = −4.47 eV) to form structure 26, which is phenol. Alternatively, cyclohexadienone can be deoxygenated directly via R50 to obtain structure 25 with 1.59 eV of energy. Further hydrogenation of structure 25 yields benzene (R35).

Once formed, the desorption of phenol from the surface is kinetically slow (3.49 × 10¹ s⁻¹ at 500 K) and energy demanding (E_des = 4.45 eV), which means it would be likely to accumulate on the surface. The adsorption energy, which can be calculated as the negative of the desorption energy, suggests that the adsorption of phenol on the surface is at least 2 eV stronger than that observed over transition metals.^{54,58,59,67,68} Thus, the accumulated phenol could convert further to benzene, which is considered further herein via two mechanisms. One possibility is the cleavage of the –OH group (R33) with an energy barrier of 1.50 eV, while alternatively the phenol α-position can be hydrogenated with an energy barrier of 2.06 eV to form structure 27 (R34). Kinetically, R34 is much slower than R33, with rate constants differing by ∼10³ s⁻¹. Therefore, the reaction will proceed via R33 to form structure 25, which will then be hydrogenated to form structure 29, benzene. The preference in our calculations towards dehydroxygenation of phenol rather than hydrogenation is in strong agreement with observations involving noble metal catalysts,^28,53,69 and in non-catalytic works.⁷⁰

The several smaller radicals formed during the outlined reaction mechanism have been considered, as they will react together to form stable entities and desorb from the surface. In our work, R37 and R38 describe the formation of methane from CH₂; similarly, R40 and R41 describe the formation of methanol from OCH₂; and R39 describes the formation of water from OH and H. To comprehensively complete the catalytic cycle, R37 to R49 are used in the microkinetic modelling.

Overall, our results show that the most favorable pathway for the formation of benzene is via cleavage of the methoxy group in guaiacol (R4) to give structure 5/20, which is then hydrogenated (R12) to give phenol (Fig. 3). The hydroxyl group of the phenol is then cleaved (R33) with an energy barrier of 1.50 eV to yield structure 25, which is hydrogenated to give benzene. All the reactions become faster with increase in temperature, except the guaiacol adsorption on the surface.

Fig. 3 The reaction profile for the most favorable pathway of HDO of guaiacol over β-Mo₂C. The activation barrier for each step is given in red, in eV, and the rate constant of the conversion at 500 K is given in purple, in s⁻¹.

To further summarise our findings, the overall activation energy for the formation of different important intermediates was calculated, as shown in Table 2. The formation of benzene via a reaction process of R4–R12–R33–R35 gives an activation energy of 1.75 eV, where reaction 33 is the rate limiting step with an activation energy of 1.50 eV. In the formation of catechol via reaction R5–R15–R16, the rate controlling step is reaction R5, and the overall activation energy for the conversion is 1.02 eV. The formation of anisole via reaction R3–R11 presents an overall barrier of 1.04 eV, which is associated with the dehydroxylation of guaiacol in reaction R3. Finally, the conversion of guaiacol to phenol is most feasible via reaction R1–R9 with an overall activation barrier of 1.15 eV during the hydrogenation of the adsorbed guaiacol.

Table 2 Overall activation barrier of the formation of benzene, catechol, anisole and phenol, and the corresponding pre-exponential factor

Overall reaction	Considered path	Activation energy (eV)	Pre-exponential factor (s^−1)
Guaiacol → benzene	R4–R12–R33–R35	1.75	6.48 × 10^13
Guaiacol → catechol	R5–R15–R16	1.02	2.26 × 10^13
Guaiacol → anisole	R3–R11	1.04	4.03 × 10^13
Guaiacol → phenol	R1–R9	1.15	4.47 × 10^13

---

3.2. Microkinetic modelling

Microkinetic modelling of the reaction network described in section 3.1 has been performed, considering specifically a batch reactor model at temperatures (T) ranging from 500 K to 700 K. The rate constants obtained from the thermochemical analyses were used to write the rate equation for 192 elementary reactions (Table S4, ESI†), and the set of coupled ordinary differential equations was solved. The initial ratio of hydrogen and guaiacol was set to unity, and in large excess with respect to the number of surface sites. The reaction profile was studied for the duration of 1 s from the initialization of the reaction (i.e. t = 0 s). At the atomic level in the DFT simulations, the reactions occur very quickly due to the absence of any heat and mass transfer resistance; furthermore, since the rate constants are first order (s⁻¹), and because the magnitudes of the rate constants are much higher than the time scale considered herein, the total time considered in the study is sufficient to model the system.

As the kinetics of formation of smaller species like water, methanol and methane are much higher than the kinetics of the formation of aromatic compounds, it is observed that methane, water, methanol and vacant sites cover 99% of the catalyst surface. The aromatic species like catechol, phenol and others, on the other hand, were present in very low coverage (<1%) at all times. The formation of the small species (water, methanol and methane) and the resultant loss in the carbon balance due to their difficult condensation has also been reported in experimental works.⁷¹ Moreover, such low coverage of the species of interest during the microkinetic modelling of guaiacol is not unique²⁸ and, as the observed trends fairly representant the reactivity on the surface of the catalyst in DFT studies, they can be relied upon for insights into the surface reaction mechanism.

At 500 K, the concentration of guaiacol on the surface increases steeply at 0.02 s (Fig. 4), as the guaiacol is adsorbed on the surface, to reach 10⁻¹⁸ ML (monolayer). During the same time interval, the concentration of surface bound phenol increases to about 10¹⁰ times the concentration of the adsorbed guaiacol, suggesting that a large fraction of guaiacol immediately converts to phenol upon adsorption. The quick reaction occurs because the rate of formation of phenol is of the order 10¹² s⁻¹ in reaction 9 (Table 1), which leads to such high concentration of phenol forming in a short time interval. A corresponding increase of phenol is also observed in the gas phase, rising to 8.59 × 10⁻²⁸ Pa in 0.02 s and maintaining a steady state thereafter. Since the catechol conversion kinetics competes with phenol formation, the catechol formation rate increases, to 5.85 × 10⁻²³ ML at 0.32 s, while the phenol conversion attains a steady state (1.99 × 10⁻²⁷ ML) after 0.18 s. Hydrogenation of phenol produces benzene, which appears on the surface in a lower concentration (3.39 × 10⁻¹¹ ML) as the kinetics of reaction 33 and 34 are much slower than phenol formation. As more catechol, phenol and benzene are produced, the concentration of guaiacol declines to 4.38 × 10⁻³² ML at 0.06 s; thereafter, the concentration of guaiacol starts to decrease steadily by 0.5 × 10⁻¹ ML per 0.01 s. However, as the smaller species like methanol and methane desorb from the surface due to fast kinetics (Table S4†), a slight increase in the production of catechol on the surface is observed from 0.30 s to 0.80 s. In the same duration, the concentration of catechol in gas phase increases by 10 times due to increased surface desorption. Thereafter, all species achieve steady state.

Fig. 4 Logarithmic graphs of the concentration of guaiacol (structure 1, blue), phenol (structure 26, red), catechol (structure 13, black) and benzene (structure 29, green) (a) on the surface, and (b) in the gas phase at 500 K, 600 K and 700 K, for the reaction time of 1 s.

At 600 K, the concentration profile of all species exhibits less changes than at 500 K, due to the changes in the rate constants of all the reactions. Similar to 500 K, the concentration of phenol rises above all species at 0.02 s; however, the coverage (5.28 × 10⁻⁸ ML) is marginally lower than at 500 K. The lower guaiacol adsorption rate is because the rate of adsorption decreases with increase in temperature. As less guaiacol adsorbs to the surface, the concentration of products formed is also lower; at 0.02 s, guaiacol coverage is slightly elevated to 7.95 × 10⁻²⁸ ML, which then slowly decreases by the order of 10⁴ to achieve steady state at 0.30 s. The concentration of benzene is 4.09 × 10⁻¹⁸ ML at 0.02 s, due to quick phenol decomposition, but the concentration falls thereafter as the selectivity shifts towards catechol formation. After 0.08 s, the concentration of catechol starts to rise, coupled with a further decrease in the guaiacol concentration; the gas phase profile shows fluctuations in the benzene concentration at this time, suggesting that the rate of benzene desorption is competitive to the rate of formation of benzene on the surface. The concentration of catechol in gas phase has an increasing rate after 0.10 s, as a consequence of the increased formation on the surface. After 0.20 s, all the species start to attain steady state with the gas phase profiles of all products showing a steady marginal increase in production of value-added products after 0.5 s.

When the reaction is considered at 700 K, the concentration profile of all the compounds reach steady state very quickly. At 0.02 s, the concentration of guaiacol is expected to be much lower than at T = 500 or 600 K, due to decrease in the guaiacol adsorption rate constant (6.45 × 10⁻⁵ s⁻¹); however, the concentration is calculated as 1.26 × 10⁻²⁸ ML, which is only marginally lower than at 600 K. The higher surface coverage is because the reaction rate of almost all the reactions considered in the kinetic model increased at 700 K, leading to faster formation of products, and thus, faster product desorption from the surface. The quicker rate of reaction recreates vacant surface sites in a shorter time span than at lower T, and leads to a higher guaiacol adsorption. The higher product concentration observed in the gas phase at 0.02 s for 700 K, when compared to lower T, is the indicator of the quick desorption.

From t = 0 to t = 0.26 s, the concentration profile of guaiacol increases sharply to the maximum, and then declines steadily to 5.95 × 10⁻³⁰ ML. After 0.26 s, a steady state is maintained until 0.38 s, from which time more products are desorbed from the surface, freeing surface sites and subsequently resulting in higher phenol, catechol and benzene formation. Further analysis of the selectivity of products shows that phenol will be the most likely product, with over 99% selectivity at all times, and this is consistent at all temperatures. The concentration of other products is very low in comparison, over the surface and in the gas phase. In the literature,⁷² high selectivity of benzene (35%) followed by phenol (30%) is reported among all other aromatic/cyclic products over Mo₂C at 623 K and 27.48 bar. Blanco et al.⁷³ report over 65% selectivity for phenol at 50 bar pressure and 623 K temperature over activated carbon supported Mo, which also has strong similarities to our modelled reaction process. The high hydrogen pressure in the experimental work could be the likely reason for the reported high selectivity of the fully hydrogenated benzene over phenol. Tran et al.⁷¹ also report phenol to desorb quickly after its formation on the same surface without undergoing further conversion. As the partial pressure of hydrogen is 0.5 in our study, a higher selectivity of phenol is obtained instead of full hydrogenation to benzene.

To determine the rate controlling step for the formation of each product in the gas phase, the degree of rate control was also analysed by the method proposed by Campbell,⁷⁴ where the forward rate constants of each reaction is perturbed by 0.1% while keeping the equilibrium constant fixed, and the changes in the conversion of products are analysed. The degree of rate control for reaction i (X_RC,i) was calculated as:

where r is the rate of reaction, k_i and k_j are the forward rate constants, and K_i is the equilibrium rate constant. The conversion of guaiacol to catechol on the surface is influenced most by reaction R16, where structure 7 is converted to catechol, and by R24, where catechol is further degraded into structure 5. As the rate of all the reactions preceding the formation of catechol is much higher than R16, R16 becomes the rate limiting step with degree of rate control (DRC) of 0.72. Its subsequent desorption from the surface has the maximum DRC. The conversion to phenol is highly favourable on the surface, with no reactions significantly influencing its formation in either positive or negative way; this outcome is primarily because the formation of phenol occurs via different competing pathways, so changes in the rate constant of one reaction does not influence its formation in any significant way. The desorption of phenol, however, is the rate limiting step showing the DRC of 1.00; the same holds true for the formation of benzene, for which the DRC is 0.99 when considering desorption from the surface. The highlighted observations suggest that desorption is the rate limiting step for the aromatic species. For the smaller species like CH₄, CH₃OH and water, the reactions describing their formation in R38, R41 and R39, respectively, show the highest DRC (∼1.00).

Temperature programmed desorption (TPD). A temperature programmed desorption simulation was also conducted, where the initial guaiacol coverage on the catalyst surface was considered at 10%, 40%, 70% and 100%. The temperature of the system was increased at a rate of 10 K s⁻¹, and the concentration change of the species in the gas phase were recorded from 300 K to 700 K. No re-adsorption of the species was allowed during the analysis, in line with the experimental procedure.

The TPD model (Fig. 5) shows that the desorption of benzene occurs near 450 K at all coverages, after it is produced from phenol. Change in concentration of benzene has been reported at temperatures above 250 K for a different surface (Pt and Pt₃Sn), at similar coverage,⁷⁵ whilst benzene is known to desorb at temperatures as low as 150 K from graphene.⁷⁶ Therefore, it is likely that the carbon-liking metal terminated surface investigated here binds the benzene strongly to the surface.⁷⁷ Phenol, on the other hand, stays on the surface till 600 K and starts desorbing thereafter. The desorption of phenol is known to occur at above 450 K over other catalyst material surfaces.^78–80 At high coverage (70%), the desorption of benzene from the surface is large, as more benzene is formed on the surface from phenol. At the same coverage, catechol desorption is also observed in large proportion after 500 K, due to its formation and accumulation on the surface. Further increasing the guaiacol coverage (100%) is detrimental to reaction rates, with the desorption rate of all species reduced due to lack of sites for the generation of species on the surface.

Fig. 5 Change in gas phase concentration of guaiacol (Structure 1, blue), phenol (Structure 26, red), catechol (Structure 13, black) and benzene (Structure 29, green), shown as a plot of change in partial pressure vs. temperature, for cases when the guaiacol coverages on the surface are (a) 10%, (b) 40%, (c) 70% and (d) 100%.

4. Conclusions

The HDO of guaiacol was studied over the (100) β-Mo₂C surface. A network of reaction pathways was proposed to consider the mechanistic aspects of producing benzene from guaiacol, proceeding via phenol and catechol intermediates. The activation barriers and reaction energies were calculated using density functional theory. Through thermochemical analysis, increases in temperature were observed as being unfavorable for the surface adsorption of guaiacol, and the most favourable pathway for guaiacol HDO is initiated with demethoxylation, which is then followed by hydrogenation to yield phenol. From phenol, the most favourable pathway for benzene formation is initiated with dehydroxylation, followed by hydrogenation.

Overall, the hydrogenation of molecules like guaiacol and phenol were observed to have higher energy barriers than the cleavage of functional groups. The carbon atoms of the functional group of the molecule, which were sp² or sp hybridized, were also noted as dissociating spontaneously, with minimum energy barriers, and the dissociated moieties occupy carbon deficient sites on the surface.

Complementing the thermochemical analysis, microkinetic modelling of the system shows that although the adsorption of guaiacol decreases with increasing temperature, higher temperatures are more favorable for further upgrading of guaiacol into valuable products such as phenol and benzene. Simulation of temperature programmed desorption shows that the optimum surface coverage of guaiacol is ∼70%, as the change in the concentration of all the products was highest in the gas phase.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thanks Igor Kowalec for technical discussions. KA is a Commonwealth Scholar, funded by the UK government. The authors thank the U. K. High Performance Computing “Materials Chemistry Consortium” (EP/R029431/1) for access to the ARCHER National Supercomputing Service, Supercomputing Wales for access to the Hawk HPC facility, which is part-funded by the European Regional Development Fund via the Welsh Government, and GW4 and the UK Met Office for access to the Isambard UK National Tier-2 HPC Service (EP/P020224/1). AJL acknowledges funding by the UKRI Future Leaders Fellowship program (MR/T018372/1).

References

1. A. V. Bridgwater and G. V. C. Peacocke, Renewable Sustainable Energy Rev., 2000, 4, 1–73.
2. A. R. K. Gollakota, M. Reddy, M. D. Subramanyam and N. Kishore, Renewable Sustainable Energy Rev., 2016, 58, 1543–1568.
3. H. Wang, J. Male and Y. Wang, ACS Catal., 2013, 3, 1047–1070.
4. G. W. Huber, I. Sara and A. Corma, Chem. Rev., 2006, 2, 4044–4098.
5. Y.-C. Lin and G. W. Huber, Energy Environ. Sci., 2009, 2, 68–80.
6. P. M. Mortensen, J. D. Grunwaldt, P. A. Jensen, K. G. Knudsen and A. D. Jensen, Appl. Catal., A, 2011, 407, 1–19.
7. J. Han, J. Duan, P. Chen, H. Lou and X. Zheng, Adv. Synth. Catal., 2011, 353, 2577–2583.
8. H. Xin, A. Vojvodic, J. Voss, J. K. Nørskov and F. Abild-Pedersen, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 115114.
9. J. K. Nørskov, F. Abild-Pedersen, F. Studt and T. Bligaard, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 937–943.
10. E. Toyoda, R. Jinnouchi, T. Hatanaka, Y. Morimoto, K. Mitsuhara, A. Visikovskiy and Y. Kido, J. Phys. Chem. C, 2011, 115, 21236–21240.
11. Y. Zhou, Z. Zhou, R. Shen, R. Ma, Q. Liu, G. Cao and J. Wang, Energy Storage Mater., 2018, 13, 189–198.
12. F. Ando, T. Tanabe, T. Gunji, S. Kaneko, T. Takeda, T. Ohsaka and F. Matsumoto, ACS Appl. Nano Mater., 2018, 1, 2844–2850.
13. H. Tao, S. Liu, J.-L. Luo, P. Choi, Q. Liu and Z. Xu, J. Mater. Chem. A, 2018, 6, 9650–9656.
14. J. Zhou, W. An, Z. Wang and X. Jia, Catal. Sci. Technol., 2019, 9, 4314–4326.
15. Z. Chen, Y. Song, J. Cai, X. Zheng, D. Han, Y. Wu, Y. Zang, S. Niu, Y. Liu, J. Zhu, X. Liu and G. Wang, Angew. Chem., Int. Ed., 2018, 57, 5076–5080.
16. J. A. van Bokhoven and J. T. Miller, J. Phys. Chem. C, 2007, 111, 9245–9249.
17. F. Morteo-Flores, J. Engel and A. Roldan, Philos. Trans. R. Soc., A, 2020, 378, 20200056.
18. R. Moreira, E. Ochoa, J. L. Pinilla, A. Portugal and I. Suelves, Catalysts, 2018, 8(4), 127.
19. R. N. Olcese, M. Bettahar, D. Petitjean, B. Malaman, F. Giovanella and A. Dufour, Appl. Catal., B, 2012, 115–116, 63–73.
20. T. Nimmanwudipong, R. C. Runnebaum, D. E. Block and B. C. Gates, Energy Fuels, 2011, 25, 3417–3427.
21. J. Sun, A. M. Karim, H. Zhang, L. Kovarik, X. S. Li, A. J. Hensley, J.-S. McEwen and Y. Wang, J. Catal., 2013, 306, 47–57.
22. D. C. Elliott and T. R. Hart, Energy Fuels, 2009, 23, 631–637.
23. C. Zhao, J. He, A. A. Lemonidou, X. Li and J. A. Lercher, J. Catal., 2011, 280, 8–16.
24. A. Gutierrez, R. K. Kaila, M. L. Honkela, R. Slioor and A. O. I. Krause, Catal. Today, 2009, 147, 239–246.
25. Y.-C. Lin, C.-L. Li, H.-P. Wan, H.-T. Lee and C.-F. Liu, Energy Fuels, 2011, 25, 890–896.
26. V. N. Bui, G. Toussaint, D. Laurenti, C. Mirodatos and C. Geantet, Catal. Today, 2009, 143, 172–178.
27. M. A. González-Borja and D. E. Resasco, Energy Fuels, 2011, 25, 4155–4162.
28. J. Lu, S. Behtash, O. Mamun and A. Heyden, ACS Catal., 2015, 5, 2423–2435.
29. R. Ma, K. Cui, L. Yang, X. Ma and Y. Li, Chem. Commun., 2015, 51, 10299–10301.
30. A. L. Jongerius, R. W. Gosselink, J. Dijkstra, J. H. Bitter, P. C. A. Bruijnincx and B. M. Weckhuysen, ChemCatChem, 2013, 5, 2964–2972.
31. M. Siaj, C. Reed, S. T. Oyama, S. L. Scott and P. H. McBreen, J. Am. Chem. Soc., 2004, 126, 9514–9515.
32. B. Martínez, F. Viñes, P. H. McBreen and F. Illas, J. Catal., 2021, 393, 381–389.
33. T. Epicier, J. Dubois, C. Esnouf, G. Fantozzi and P. Convert, Acta Metall., 1988, 36, 1903–1921.
34. T. Wang, X. Tian, Y. Yang, Y.-W. Li, J. Wang, M. Beller and H. Jiao, Surf. Sci., 2016, 651, 195–202.
35. K. Xiong, W. Yu, D. G. Vlachos and J. G. Chen, ChemCatChem, 2015, 7, 1402–1421.
36. L. M. Falicov and G. A. Somorjai, Proc. Natl. Acad. Sci. U. S. A., 1985, 82, 2207–2211.
37. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun., 2009, 180, 2175–2196.
38. A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. Bjerre Jensen, J. Kermode, J. R. Kitchin, E. Leonhard Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
39. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
40. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
41. M. J. Frisch, J. A. Pople and J. S. Binkley, J. Chem. Phys., 1984, 80, 3265–3269.
42. K. Agrawal, A. Roldan, N. Kishore and A. J. Logsdail, Catal. Today, 2022, 384, 197–208.
43. J. Nocedal and S. J. Wright, in Numerical Optimization, Springer New York, 2nd edn, 2006, pp. 66–100.
44. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
45. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
46. J. A. Garrido Torres, P. C. Jennings, M. H. Hansen, J. R. Boes and T. Bligaard, Phys. Rev. Lett., 2019, 122, 156001.
47. A. Roldan, 2021, https://github.com/Roldan-Group/Kinetics.
48. X. Lu, J. Zhang, W.-K. Chen and A. Roldan, Nanoscale Adv., 2021, 3, 1624–1632.
49. I. Chorkendorff and J. W. Niemantsverdriet, Concepts of Modern Catalysis and Kinetics, 2003.
50. H. Eyring, J. Chem. Phys., 1935, 3, 107.
51. M. G. Evans and M. Polanyi, Trans. Faraday Soc., 1935, 31, 875–894.
52. K. W. Kolasinski, Surface Science, Foundations of Catalysis and Nanoscience, John Wiley & Sons, Ltd, Chichester, UK, 3rd edn, 2012.
53. J. Lu and A. Heyden, J. Catal., 2015, 321, 39–50.
54. A. M. Verma and N. Kishore, J. Mol. Model., 2018, 24, 254.
55. K. Lee, G. H. Gu, C. A. Mullen, A. A. Boateng and D. G. Vlachos, ChemSusChem, 2015, 8, 315–322.
56. H. Fang, A. Roldan, C. Tian, Y. Zheng, X. Duan, K. Chen, L. Ye, S. Leoni and Y. Yuan, J. Catal., 2019, 369, 283–295.
57. C. Chiu, A. Genest, A. Borgna and N. Rösch, ACS Catal., 2014, 4, 4178–4188.
58. S. J. Carey, W. Zhao, Z. Mao and C. T. Campbell, J. Phys. Chem. C, 2019, 123, 7627–7632.
59. J. Lee, S. Ryu and S. K. Kim, Surf. Sci., 2001, 481, 163–171.
60. J. Chang, T. Danuthai, S. Dewiyanti, C. Wang and A. Borgna, ChemCatChem, 2013, 5, 3041–3049.
61. V. N. Bui, D. Laurenti, P. Afanasiev and C. Geantet, Appl. Catal., B, 2011, 101, 239–245.
62. M. V. Bykova, D. Y. Ermakov, V. V. Kaichev, O. A. Bulavchenko, A. A. Saraev, M. Y. Lebedev and V. Yakovlev, Appl. Catal., B, 2012, 113–114, 296–307.
63. M. Nagai, H. Tominaga and S. Omi, Langmuir, 2000, 16, 10215–10220.
64. J. R. McManus and J. M. Vohs, Surf. Sci., 2014, 630, 16–21.
65. W. S. Lee, Z. Wang, R. J. Wu and A. Bhan, J. Catal., 2014, 319, 44–53.
66. Z. Cai, F. Wang, X. Zhang, R. Ahishakiye, Y. Xie and Y. Shen, Mol. Catal., 2017, 441, 28–34.
67. F. Maldonado, L. Villamagua and R. Rivera, J. Phys. Chem. C, 2019, 123, 12296–12304.
68. A. M. Verma and N. Kishore, R. Soc. Open Sci., 2017, 4, 170650.
69. M. Saleheen, A. M. Verma, O. Mamun, J. Lu and A. Heyden, Catal. Sci. Technol., 2019, 9, 6253–6273.
70. K. Agrawal, A. M. Verma and N. Kishore, ChemistrySelect, 2019, 4, 6013–6025.
71. C.-C. Tran, O. Mohan, A. Banerjee, S. H. Mushrif and S. Kaliaguine, Energy Fuels, 2020, 34, 16265–16273.
72. C. C. Tran, Y. Han, M. Garcia-Perez and S. Kaliaguine, Catal. Sci. Technol., 2019, 9, 1387–1397.
73. E. Blanco, C. Sepulveda, K. Cruces, J. L. García-Fierro, I. T. Ghampson and N. Escalona, Catal. Today, 2020, 356, 376–383.
74. C. T. Campbell, J. Catal., 2001, 204, 520–524.
75. B. M. Wong, G. Collinge, A. J. R. Hensley, Y. Wang and J.-S. McEwen, Prog. Surf. Sci., 2019, 94, 100538.
76. R. S. Smith and B. D. Kay, J. Phys. Chem. B, 2018, 122, 587–594.
77. T. Mo, J. Xu, Y. Yang and Y. Li, Catal. Today, 2016, 261, 101–115.
78. E. Castillejos-López, D. M. Nevskaia, V. Muñoz and A. Guerrero-Ruiz, Carbon, 2008, 46, 870–875.
79. B.-Q. Xu, T. Yamaguchi and K. Tanabe, Mater. Chem. Phys., 1988, 19, 291–297.
80. X. Xu and C. M. Friend, J. Phys. Chem., 1989, 93, 8072–8080.

---

Footnote

† Electronic supplementary information (ESI) available at DOI: 10.1039/d1cy01273h, including: Additional information is provided in a separate section containing details of the optimised configurations observed for all intermediates considered in the simulations, and the data derived from the microkinetic simulations. All structures calculated in this work have been uploaded to the NOMAD repository (DOI: 10.17172/NOMAD/2022.01.05-1).

---