Molecular Catalysis Density functional theory study of the zeolite-mediated tautomerization of phenol and catechol

Because the structure of lignin consists mostly of inter-linked phenolic monomers, its conversion into more valuable chemicals may beneﬁt from isomerization processes that alter the electronic structure of the aromatic rings. The tautomerization of phenolic-type compounds changes the hybridization from sp 2 to sp 3 of the carbon atom at the ortho position, which disables the aromaticity and facilitates the subsequent hydrogenation process. Here, we have performed a Density Functional Theory study of the tautomerization of phenol and catechol at the external surface of zeolite MFI. The tautomerization starts with the adsorption of the molecule on three-coordinated Lewis acid sites, followed by the dissociation of the phenolic hydroxyl group, with the transfer of the proton to the zeolite framework. The rotation of the deprotonated molecule enables a more favourable orientation for the back-transfer of the proton to the carbon atom at the ortho position. The energy barriers of the process are smaller than 55 kJ/mol, suggesting that this transformation is easily accessible under standard reaction conditions.


Introduction
Lignin is one of the most abundant components of biomass, representing 10-35% in mass and 40% in energy [1,2]. However, lignin is a three-dimensional polymer with a high number of C C and C O linkages between the phenolic building blocks, which limits its use as a source of more valuable chemicals and fuels [2][3][4]. Various processes have been reported to depolymerize lignin to smaller and soluble components, which are then chemically transformed depending on their future applications. Relevant examples include the mechano-catalytic depolymerisation of lignocellulose [5], the catalytic fast pyrolysis [6] or the metal-catalysed hydroprocessing of organosolv lignin [7]. Lignin-derived compounds have significant potential to replace fossil fuels as a source of heat and power [2], but first their energy content has to be enhanced by increasing the C:O and H:C ratios [1,8].
The modification of the C:O and H:C ratios is primarily restricted by the aromatic moieties that constitute the lignin: the C O bond in phenolic species is approximately 84 kJ/mol stronger than the same type of bond in aliphatic compounds [9]. In this regard, mono-and bi-metallic nanoparticles supported on zeolites show promising results towards the hydrodeoxygenation (HDO) of lignin [7,[10][11][12][13]. These methods effectively increase the C:O and H:C ratios whilst using milder conditions than pyrolytic processes [6,[14][15][16], leading to higher yields and selectivity with practically no solid residues. The role of the zeolite in the HDO treatment is assumed to be related to the removal/addition of aliphatic side chains and the dehydration of aliphatic alcohols, which is carried out by Brønsted and Lewis acid sites. The hydrogenation/hydrogenolysis of the aromatic ring is essentially considered to be due to the metal catalyst, although an enhancement of its function is noted in the presence of the zeolite [13].
In the present work, we have analysed the activation of phenolic monomers on the external Lewis acid sites of zeolite MFI. The proposed transformation is the keto-enol tautomerization shown in Scheme 1, which involves the transfer of the hydroxyl hydrogen of phenol to one of the carbon atoms at the ortho position, thereby changing the hybridization of the carbon atom from sp 2 to sp 3 . This loss of molecular aromaticity promotes the hydrogenation of the molecule [17,18]. For instance, the initial tautomerization of m-cresol to an unstable cyclo-hexadienone has been reported to be the preferred pathway towards HDO over the Pt/SiO 2 catalyst [19]. Likewise, the tautomerization of phenol is considered a possible intermediate process during the hydroprocessing on Pd/SiO 2 , Pd/Al 2 O 3 , and Pd/ZrO 2 [20]. Similar bifunctional catalysts use zeolite as a support during the HDO of lignin and phenols [10,13], which has prompted us to use computational tools to understand the potential role of the zeolite's active sites in the tautomerization process.
Previous reports have underlined the synergy arising at the metal-support interface for hydrogenation processes, derived from the combined action of support and nanoparticles; this cooperation rules the pathway and selectivity of the reaction during the hydroprocessing [20,21]. Metal nanoparticles should be found at the external surface of the zeolite because their sizes are larger than 2 nm [10,13], which is bigger than the pore dimension of zeolite frameworks. As such, the external active sites of zeolites will be relevant in the synergistic effect of the metal-support interface. Accordingly, we have focused our attention in this study on the adsorption and reaction of representative molecules, phenol and catechol, on the external surface of the zeolite MFI, modelled here as a periodic slab. This model is not far from reality; MFI nanosheets have been successfully synthetized as thin as a single unit cell along the [010] direction, whilst extending two-dimensionally over its (010) plane [22]. We have chosen the three-coordinated aluminium centre within the zeolite framework as the active site to adsorb the molecules. This site presents Lewis acidity owing to the undercoordination of the aluminium atom, as shown in Scheme 2, which has been shown by computer simulations to be thermodynamically and kinetically favoured [23,24].
The first part of the result section is dedicated to the determination of the most stable aluminium-substituted silanol to adsorb the molecules. In the second part, we have presented a detailed mechanistic analysis of the adsorption and reaction of phenol, including the energy barriers along the pathway, to achieve its isomerization. We have examined the catechol adsorption in the final section, and analysed the consequences of the second hydroxyl group in the tautomerization process.

Computational methods
We have performed all calculations within the framework of the Density Functional Theory, using the plane wave code Vienna Ab-initio Simulation Package (VASP) [25][26][27][28]. We used the generalized gradient approximation (GGA), as derived by Perdew, Burke and Ernzerhof (PBE) [29], in combination with Grimme's correction for the dispersion forces to account for the exchange-correlation energy at short-and long-range. In a previous report, we have shown that the inclusion of Grimme's atomic pairwise dispersion correction (PBE-D2) [30] decreases the error of the computed volume of the MFI unit cell, and also gives a better prediction of the bulk modulus [24]. In the present work, we have considered a more recent approach to the dispersion interactions, which incorporates geometry information into the ab-initio parametrization of the dispersion coefficients (PBE-D3) [31,32]. The PBE-D3 method closely agreed with PBE-D2, but the best outcome was achieved considering the Becke-Johnson damping (PBE-D3 BJ) protocol, reducing the volume error of the unit cell to only 0.5% (see Table 1). This result prompted us to choose the PBE-D3 BJ approach to perform all calculations.
We have optimized the orthorhombic unit cell of the MFI framework, available from the Structure Database of the International Zeolite Association (IZA) [35]. A representation of the optimized MFI framework is shown in Fig. 1. A set of fixed-volume calculations were performed, allowing the relaxation of the lattice shape and atomic positions. Thereafter, the correlation of the optimized lattice energy versus the cell volume was fitted to the Birch-Murnaghan equation of state [34]. This methodology eliminates the problems related to the Pulay stress [36] whilst giving the equilibrium volume and the bulk modulus as adjustable parameters (see Table 1). Table 1 Structural parameters after optimization of the bulk unit cell of the zeolite MFI with the PBE functional in combination with the Grimme's methods (D2 and D3) that account for the Van der Waals interactions. The valence electrons were explicitly treated by means of a plane wave basis set with energy cutoff of 550 eV; their nodal features and interaction with the inner core were described by the projectedaugmented-wave method (PAW) [37,38]. The Gaussian smearing method was included to improve the electronic convergence, setting band widths of 0.1 eV (zeolite and molecule/zeolite) or 0.01 eV (isolated molecules) [39,40]. Owing to the large dimension of the MFI unit cell, with lattice vectors of over 13 Å, only the Gamma point was used to perform the numerical integration within the Brillouin zone. The convergence criteria for the electronic and ionic forces were set at 10 −5 eV and 0.03 eV/Å, respectively.
The binding free energies of the molecules to the acid sites were calculated according to the following equation: where F molecule/zeolite is the free energy of a system constituted by the zeolite and a molecule of phenol/catechol adsorbed on the Lewis acid site, F zeolite is the free energy of the bare zeolite, and F gas molecule is the energy of a molecule of phenol/catechol in the gas phase. The term gas phase is used for a single molecule optimized within a 30 × 31 × 32 Å 3 supercell. The values of the free energies were calculated after including the zero-point energy correction and the entropic effects at 473 K, which is a fair representation of the range of temperatures commonly used [10][11][12]41]. Only the vibrational contribution to the entropy was considered for the zeolite/adsorbate systems. In that case, the harmonic vibrational partition function was calculated using the normal modes derived from the movement of the Al atom, the first sphere of SiO 4 tetrahedra binding the Al and the adsorbed molecule of phenol or catechol, whilst the rest of the atoms were kept frozen. In the case of molecules in the gas phase, the rotational and translational contributions to the entropy were also considered, assuming them to be ideal gases.
We applied the Nudged Elastic Band (NEB) method to each reaction step to determine the initial guess of the transition state (TS) [42,43]. This technique starts the exploration of the minimum energy pathway (MEP) between reactants and products by linking both states with evenly spaced images connected by harmonic springs. The perpendicular component of the spring force, which tends to make the images collinear, is zeroed. The same treatment is given to the parallel component of the potential energy gradient, which tends to move the images towards the closest minimum. In this work, we have used between one and four images to create the elastic bands, depending on the particular configuration, and a value of 0.1 eV/Å was used as convergence criterion for the adjusted force on each image. The position of the TS was refined using the improved dimer method (IDM) [44], which increases the robustness and efficiency of the original dimer method (DM) proposed by Henkelman and Jónnson [45]. The DM is an algorithm to search for local saddle-points using the potential energy gradients, thus avoiding expensive evaluations of the Hessian matrix during the search. Two configurations in the multi-dimensional space, separated by a small vector (dimer axis), are used to follow the mode with the lowest negative curvature in the potential energy. The starting path to align the dimer is the hardest imaginary mode of the TS guess derived from the NEB calculation. In this way a faster convergence is achieved. The IDM is considered converged when a gradient criterion is reached along the lowest negative curvature. In the present work, we have used a threshold of 0.03 eV/Å. During the geometry optimization of the TS, using the NEB method and IDM, only the Al atom, the first sphere of the SiO 4 tetrahedra binding the Al, including any H binding the framework O atoms, and the molecule were allowed to relax, whereas the rest of the atoms were kept frozen in their positions. The systems containing the reactant and product of each step were fully relaxed without any geometry constraints.
The images of the structural geometries were obtained with the code Visualization for Electronic and Structural Analysis (VESTA 3) [46].

Slab model and aluminium distribution
We used the optimized MFI unit cell to construct the zeolite slabs by keeping the periodicity along the [100] and [001] directions, and cutting along the [010] direction, which coincides with the straight pore orientation. Silanol groups were formed after saturation of the dangling Si-O bonds at the external surface. A vacuum layer of 20 Å was inserted to avoid interactions between periodic images. Thereafter, we performed the Al substitution at the twelve non-equivalent T-sites exposed at the external surface. Two slabs were tested: the first one was constructed with two pentasil layers, which reproduced the thinnest experimental thickness of an MFI nanosheet ( Fig. 2a) [47], whilst the second model was formed by a single pentasil layer in order to reduce the computational cost of the calculations (Fig. 2b). Both thicknesses were compared during the analysis of the Al substitution to validate the use of the one-pentasil slab model.
The Al substitution creates a negative charge within the framework, owing to the replacement of Si atoms with formal charge 4+ by Al with 3+. This negative charge is balanced by tetraalkylammonium cations, which also act as structure directing agents (SDA) during the synthesis of Al-doped zeolite [48]. We considered the explicit presence of the SDA by placing a tetramethylammonium (TMA) cation in the centre of the interception between the sinusoidal and straight channels, as shown in Fig. 3, where we have disregarded the fact that in experiment the alkyl chain is much longer than one carbon atom [48]. However, we decided to use the tetramethylammonium to simplify the explicit consideration of the counter-ion.
Despite the difference in thickness, the two slabs shown in Fig. 2 led to practically the same order of stability among the Alsubstituted T-sites, as shown in Fig. 4a. Furthermore, the distortion of the straight pore followed a similar pattern of expansions and contractions for the twelve different Al substitutions, regardless of the size of the slab (see Fig. 4b and c). Having ascertained that the one-pentasil slab adequately reproduced the results of the two- pentasil slab, we have used the former to analyse the adsorption and tautomerization of phenol and catechol.
The determination of the Al distribution is a challenging task, both experimentally and theoretically, due to the large number of factors influencing the synthesis process of the final zeolite [49][50][51]. However, the more stable structures will occur more frequently within the zeolite framework [51]. According to the results of the Al distribution in the one-pentasil slab (see Fig. 4a), the T9 site was the most stable Al-substituted silanol out of four available positions: T7, T9, T10 and T12 (see Fig. 2c). The dehydration of Al-substituted silanols has an energy barrier of between 10 and 20 kJ/mol, leading to a water molecule adsorbed on a three-coordinated Al centre, which is more stable than the initial Al-substituted silanol by 12-45 kJ/mol [23,24]. The water that is formed during dehydration may be removed after thermal pre-treatment, leaving a naked three-coordinated Lewis site (see Fig. 5) [52]. We have therefore used the dehydrated form of the T9 site to analyse the adsorption of phenol and catechol on a Lewis site at the external surface of zeolite MFI. Because the Al-substituted T9 site is three-coordinated after assuming its dehydration, we did not have to consider any longer the presence of the counter-ion (TMA) within the system.

Adsorption of phenol
We have considered two different configurations for the adsorption of phenol on the Al-substituted T9 site, referred to as non-planar (Fig. 5a) and co-planar (Fig. 5b) configurations. We have also examined the binding through the C atoms at ortho and para positions (C2 and C4 in Fig. 5). Fig. 5 shows the most stable geometries, after optimization, for the adsorption of phenol through its O atom with non-and coplanar configurations (hereafter we use O ph and H ph to refer to the  hydroxyl group of phenol). The co-planar adsorption of phenol on the Lewis site T9 was stronger than the non-planar configuration, with binding free energies of −42 and −13 kJ/mol, respectively. The binding free energy of phenol through the C2 and C4 atoms on T9 was +4 and −7 kJ/mol, respectively. Hence, the adsorption of phenol takes place preferentially through the O ph atom, and adopting a co-planar configuration.
The elongation of the O ph C bond was the main deformation of phenol upon adsorption, becoming longer for stronger O ph -Al interactions. For instance, the average O ph C bond length reached values of 1.407 and 1.431 Å for the non-and co-planar configurations, respectively, compared with a gas phase value of 1.376 Å.

Dissociation of phenol
We have considered the dissociation of phenol and the formation of phenoxide as the first elemental step towards the tautomerization (see Scheme 1). In this process, the O ph H ph bond was broken and H ph transferred to the closest framework O atom binding the Al. We have used the geometries with the strongest binding free energies on T9 to investigate the formation of phenoxide. The Bader analysis of atomic charges [53][54][55] confirmed that   6. shows the optimized structures with non-and co-planar configurations after the dissociation of the O ph H ph bond at the T9 Lewis acid site. The relative energy of the system increased by 11 and 18 kJ/mol, with the O ph C bond length decreasing to 1.354 and 1.369 Å upon deprotonation, for the non-and co-planar configurations, respectively (see Fig. 6.).
The energy barrier for the formation of phenoxide on T9 is practically independent of the phenol configuration. We obtained values of 41 and 49 kJ/mol, respectively, depending on whether the molecule had a non-planar or a co-planar configuration (see Fig. 7.). These barriers are at least 20 kJ/mol higher than the dehydration process of an Al-substituted silanol, which involves the deprotonation of the framework O atom binding the Al (Brønsted acid) and the transfer of the proton to the hydroxyl group attached to the Al, with the subsequent formation of water (see Scheme 2) [23,24]. In comparison, an activation energy of 19 kJ/mol has been reported for the dissociation of phenol on Pd(111), increasing to 46 kJ/mol when the adsorption takes place on Pt(111) [56]. In addition, the dissociation of the O H bond of m-cresol on Pt(111) and Ru(0001) reveals barriers of 39 and 45 kJ/mol, respectively [57]. Therefore, our calculated activation energies show that the dissociation of phenol on the Lewis acid sites of zeolites is as probable as on metallic surfaces, and thus can take place at relatively low temperatures.
The reverse reaction, which regenerates phenol from phenoxide, had barriers of approximately 30 kJ/mol. However, the phenoxide is free to rotate, and change its orientation, whilst H ph remains bound to the framework O atom. This alternative route is energetically less expensive than the regeneration of phenol, and has the additional advantage of favouring the transfer of H ph to C2.

Rotation of dissociated phenol
The isomerization of phenol to form cyclohexa-2,4-dien-1-one (see Scheme 1) requires the H ph transfer to one of the C atoms at the ortho position. In the initial relaxed orientation of phenoxide adsorbed on the Lewis acid T9 site (Fig. 6.), the direct transfer of the proton from O8 to C2, forming cyclohexa-2,4-dien-1-one, was less probable than the regeneration of phenol. However, after a change in the phenoxide orientation through the rotation of the dihedral angle C2-O ph -Al-O, shown in Fig. 8, the transfer of H ph became feasible. Fig. 9 shows the optimized configurations after the phenoxide rotation. The non-planar adsorption had a binding free energy of −22 kJ/mol, whilst the co-planar adsorption was −27 kJ/mol (see Fig. 9). The new orientations led to a decrease in the H ph C2 distances to approximately 2.0 Å and a small increase in the H ph O8 bond length from 0.99 Å to approximately 1.01 Å. The elongation of the H ph O8 bond indicates the interaction between H ph and the aromatic states of phenol, which favours the proton transfer to form cyclohexa-2,4-dien-1-one (see Fig. 9).
The calculated energy barrier of the rotation pathway was 12 kJ/mol for the co-planar configuration. In the case of the nonplanar configuration, the reorientation of the phenoxide occurred smoothly, without measurable barrier, as shown in Fig. 7. The newly adopted orientations of the aromatic ring hinder the regeneration of phenol, for which a reverse rotation is needed, but our observations indicate that phenoxide is likely to rotate constantly around the axis O ph -Al.

Formation of cyclohexa-2,4-dien-1-one
The final step of the phenol tautomerization involves the transfer of H ph to C2 with the subsequent formation of cyclohexa-2,4-dien-1-one. The optimized structures after the proton transfer are shown in Fig. 10. The non-planar configuration of the newly formed cyclo-hexadienone was more stable than the co-planar orientation by 8 kJ/mol. The non-planar configuration had a binding free energy of −27 kJ/mol, with reference to a molecule of phenol in the gas phase together with the bare zeolite. The cyclo-hexadienone in the gas phase is less stable than phenol by 70 kJ/mol, which means that once the tautomerization has occurred, regeneration of phenol is still more likely than desorption of cyclo-hexadienone from the Lewis acid site, although solvent effects may encourage the desorption process. However, the surface-bound cyclohexa-2,4-dien-1-one may undergo further hydrogenation without regenerating phenoxide or phenol [20,21].
After the formation of cyclohexa-2,4-dien-1-one, the O ph -Al bond length increased from approximately 1.75 Å in the phenoxide to 1.85 Å, which was accompanied by a decrease of the O ph C bond length to an average of 1.29 Å in both the non-and co-planar configurations.
We have calculated an energy barrier for the transfer of H ph to C2 of 23 and 27 kJ/mol for the two configurations (see Fig. 7.), which is almost half the barrier for the phenoxide formation. This barrier is significantly lower than the equivalent process for phenol and m-cresol on Pt(111), Pd(111) and Ru(0001), where values between 63 and 93 kJ/mol have been reported [56,57]. At the same time, an activation energy of 37 kJ/mol has been calculated for the waterassisted tautomerization of phenol at the Ru-TiO 2 (110) interphase [21], which is within our proposed range of activation energies. We also evaluated an energy barrier of 242 kJ/mol for the direct tautomerization of phenol in the gas phase, which was more than four-fold larger than the barriers along the pathways depicted in Fig. 11. Configurations of the catechol molecule. Configuration (a) is more stable than (b) by 67 kJ/mol (including zero-point energy and entropy corrections at 437 K). Fig. 7. The magnitude of the activation energy in the gas phase agreed with previous theoretical works that portray the energetics of this process as highly demanding [58,59]. Our results show that the tautomerization of phenol is effectively catalysed on the external Lewis acid sites of zeolite MFI, becoming kinetically viable under relative low temperature conditions.

Adsorption of catechol
In the gas phase, catechol had two orientations for the OH groups, as shown in Fig. 11. The most stable geometry (by 67 kJ/mol), where one OH group is oriented towards the second (shown in Fig. 11a), was adsorbed on the T9 site following the same protocol as for phenol, i.e., placing one of the two O ph close to the Al atom; this arrangement is referred along the text as the monodentate adsorption (MOA). We used the nominative second O 2ph to refer to the O ph not binding the Al centre, and not going through dissociation during the first step of tautomerization. We have also utilized the less stable configuration of catechol, with the two O ph H ph bonds oriented in opposite directions (see Fig. 11b), to simultaneously bind both O ph atoms to the Al centre; this arrangement is referred to as the bidentate adsorption (BIA). In addition, both non-and co-planar configurations were considered for each adsorption.
The structures with the strongest binding free energies within each configuration are shown in Fig. 12, (the energy of the most stable structure of catechol in the gas phase was used as a reference for both the MOA and the BIA). In the case of the MOA, binding free energies of −17 and −40 kJ/mol were calculated for the non-and co-planar configurations, respectively (see Fig. 12a and b). The O 1ph -Al distance decreased to ∼2.0 Å during relaxation, whilst the second O 2ph formed an H-bond with a nearby framework O atom, with O 2ph H 2ph · · ·O distance of ∼1.7 Å. At the same time, the BIA with non-planar configuration maintained both O ph atoms bound to the Al centre during relaxation at the expense of an energy penalty that weakened the adsorption, rendering a binding free energy of +6 kJ/mol (see Fig. 12c). However, the BIA with co-planar configuration relaxed to an MOA configuration (see Fig. 12d), with a binding free energy of −39 kJ/mol (for consistency, we labelled this structure as BIA with co-planar configuration, although knowing that it relaxed to a MOA after optimization). Therefore, the adsorption strength and orientation of catechol seems to follow the same trends of phenol, i.e., through a monodentate adsorption with co-planar configuration, whilst the second O 2ph forms stabilising H-bonds with nearby framework O atoms. However, the relative relevance of the BIA of catechol along the tautomerization process increases by kinetic and thermodynamic factors, as discussed below. Table 2 shows the binding free energy of the species involved in the tautomerization process, starting from the most stable geometries shown in Fig. 12. In contrast with phenol, the energy of the systems catechol/T9 decreased by 4-29 kJ/mol after the dissociation of the O 1ph H 1ph bond, which was expressed by binding free energies that ranged from −26 to −45 kJ/mol (see fourth column of Table 2 and Fig. 13).

Dissociation of catechol
The BIA took precedence over the MOA after the dissociation of the O 1ph H 1ph bond in response to the elongation of the Al-O8 bond, as shown in Fig. 13. The tetrahedron was transformed into a trigonal bipyramid, with the base formed by the dissociated O 1ph and the two non-protonated framework O atoms binding the Al, whilst the second O 2ph and the framework O atom binding H 1ph were placed at the top of the pyramids (see Fig. 13). The MOA with co-planar configuration was the only structure that did not adopt a BIA after dissociation, owing to an unfavourable orientation that hindered the interaction between the second O 2ph and the Al (see Fig. 13b). Therefore, although the BIA was not energetically favourable before the O 1ph H 1ph bond dissociation, it became Table 2 Binding free energies (in kJ/mol) of each of the species involved in the tautomerization of catechol starting from the structures in Fig. 12 (the barrier for the ring rotation is not included). a Labelling for the adsorbed geometries. The first three letters specify if catechol had a pre-optimized geometry that followed a MOA or a BIA, the fourth and fifth letters specify if the molecule was adsorbed with non-planar (NP) or co-planar (CP) configuration. The sixth and seventh letters specify which C atom bound the proton after the tautomerization: the one binding an H atom (CH) or the one binding the second O 2ph (CO).
b System after catechol adsorption with the strongest binding free energy. c Transition state for the dissociation of the O 1ph H 1ph bond. The activation energy is given within parenthesis in kJ/mol. d System after the dissociation of the O 1ph H 1ph bond. e System after the rotation of the dissociated catechol molecule. f Transition state for the formation of the C H 1ph bond. The activation energy is given within parenthesis in kJ/mol. g Catechol tautomer. h The structure relaxed as a co-planar MOA during the optimization of non-dissociated phenol (see Fig. 12d).
likely after the H 1ph transfer, allowing the simultaneous interaction of the two O ph of catechol with the Al centre. The dissociation of the O 1ph H 1ph bond in catechol tended to have smaller activation energies than phenol, by at least 5 kJ/mol, when the values in the third column of Table 2 are compared with those in Fig. 7. The lower energy barriers were due to the stabilising effect of the second O 2ph atom. During the dissociation, the Al-centred tetrahedron was distorted as a consequence of the shortening of the distance between the O 1ph atom and the framework O8 atom, reaching a minimum in TS 1 (see Fig. 14) and thus decreasing the distance that H 1ph needed to span to move from catechol to the zeolite. The second O 2ph stabilised the structure by getting closer to the Al atom during the O 1ph H 1ph dissociation in order to compensate for the distortion of the Al-centred tetrahedron, and eventually finished with the adoption of a BIA structure (see Fig. 13). For example, in the TS 1 of the MOA with non-planar configuration (see Fig. 14a), the H-bond O 2ph H 2ph · · ·O18 was broken, and the O 2ph H 2ph bond reoriented to allow the interaction between O 2ph and Al, decreasing the O 2ph -Al distance from 3.255 Å to 2.261 Å. This reorientation lowered the energy of the TS 1 , with both O ph remaining bound to the Al, leading to a dissociation barrier for catechol of 34 kJ/mol. An exceptionally low activation energy of 12 kJ/mol was calculated for the BIA with non-planar configuration (BIA-NP-CH/CO entries in third column of Table 2), as a consequence of the similarity between the structures of the pre-dissociated catechol (Fig. 12c) and the TS 1 (Fig. 14c), with both configurations having the two O ph bound to the Al centre. The MOA system with co-planar configuration (shown in Fig. 14b) was the only structure where the stabilising effect of the second O 2ph could not take place, owing to an unfavourable molecular orientation that hampered the interaction between the second O 2ph and the Al centre (see Fig. 13c); as a consequence, the energy barrier increased to 55 kJ/mol, which is at least 6 kJ/mol higher than in phenol.
After analysing the importance of the BIA in the TS 1 and product of the O 1ph H 1ph dissociation, the preferred configuration for the catechol adsorption should be reconsidered. We stated above that the MOA of catechol from the gas phase is more stable than BIA by at least 23 kJ/mol (see Fig. 12 and second column of Table 2). However, catechol may adsorb from the gas phase through a less thermodynamically favourable BIA with non-planar configuration, to rapidly dissociate by overcoming a low energy barrier (see BIA-NP-CH/CO entries in Table 2), and conserving the bidentate adsorption at the end of the dissociation.

Formation of catechol-derived cyclo-hexadienones
In contrast with phenol, the C atoms at ortho positions in catechol are not equivalent: one binds the second O 2ph (referred as C2) and the other binds an H atom (referred as C6). The transfer of H 1ph to C2 generated cyclohexa-3,5-dien-2-hydroxy-1-one (shown in Fig. 15a), which had unfavourable binding free energies that ranged between +21 and +47 kJ/mol (see seventh column of Table 2). On the other hand, when H 1ph was transferred to C6, with the subsequent formation of cyclohexa-2,4-dien-2-hydroxy-1-one (shown in Fig. 15b), a more thermodynamically stable product was obtained, with binding free energies between −23 and −12 kJ/mol (see seventh column of Table 2), ruling in this way the selectivity of catechol tautomerization.
The thermodynamic preference for the formation of cyclohexa-2,4-dien-2-hydroxy-1-one (Fig. 15b) over cyclohexa-3,5-dien-2hydroxy-1-one (Fig. 15a) is explained by considering the conjugation between the electronic lone pair of the second O 2ph and the conjugated double bonds that remain in the carbon ring, as shown in Scheme 3. When H 1ph is transferred to C2, the conjugation, which increases the stability of the molecule, is disrupted. In contrast, the H 1ph transfer to C6 leaves the conjugation untouched together with its stabilising effect. Furthermore, as a consequence of the formation of the cyclo-hexadienone, the C atom at the ortho position changes its hybridization from sp 2 to sp 3 , which moves the group (an aromatic H or the second O 2ph ) connected to the C atom out of the ring plane. However, in the case of the H ph transfer to C2, the second O 2ph tended to preserve its position on the ring plane, causing structural distortions. The reluctance of the second O 2ph to abandon its initial orientation forced the C2-H 1ph bond to be almost perpendicular to the ring plane, with a dihedral angle H 1ph C2-C1-O 1ph of 81 • , as shown in Fig. 15a. In comparison, the angle H 1ph C6-C1-O 1ph maintained a favourable value of 56 • (see Fig. 15b), which contributed to the higher stability of cyclohexa-2,4-dien-2hydroxy-1-one against cyclohexa-3,5-dien-2-hydroxy-1-one. The same trend was followed by the TS 2 of each alternative pathway: the energy barriers towards cyclohexa-3,5-dien-2-hydroxy-1-one were at least 30 kJ/mol larger than the equivalent values for the formation of cyclohexa-2,4-dien-2-hydroxy-1-one (see sixth column of Table 2). A similar explanation is provided for this outcome: H 1ph started to disrupt the conjugation of O 2ph at the TS 2 .
We have shown that the tautomerization of phenol and catechol is effectively catalysed by the three-coordinated Lewis acid site T9 at the external surface of zeolite MFI. This process is easily attainable energetically and involves only proton transfers and the rotation of the aromatic ring. The Lewis acidity is essential for this transformation, facilitating the adsorption of the molecule, the dissociation of the O ph H ph bond and the formation of the C H bond that culminates in the tautomerization. This suggests that metal substitutions conferring Lewis acidity to the zeolite should play an important role in the HDO of lignin-derived compounds by catalysing the tautomerization. There is a large collection of potential dopants with Lewis acid character to be tested in the HDO of lignin, ranging from the dehydrated, three-coordinated M 3+ species, such as Al 3+ , to M 4+ substitutions that do not need dehydration to behave as Lewis acids [60,61]. For instance, it has been reported that Sn-substituted zeolite BEA is an efficient catalyst for the conversion of glucose [62,63].
Several zeolites are known to have layered forms, e.g. MFI, FER and MWW [64,65], and thus have an external surface with relatively enhanced area, which can support metal nanoparticles. This two-dimensional morphology contributes to the formation of three-coordinated Lewis acids [24], and thus to the occurrence of tautomerization during the hydroprocessing. Layered zeolites have the additional advantage of making accessible a higher number of active sites to bulkier molecules that cannot enter the micropore system [66].
Finally, we cannot guarantee that the tautomerization will occur during the HDO of lignin under a comprehensive range of conditions. The effect on the tautomerization of solvents and the presence of competing transformations, as well as the effect of the tautomerization on the overall turnover of the hydroprocessing, should be further investigated, using both computational and experimental techniques.

Conclusions
We have analysed the tautomerization of phenol and catechol mediated by the Lewis acid sites at the external surface of Alsubstituted zeolite MFI by means of Density Functional Theory including dispersions corrections. Our results show that phenol and catechol adsorb with similar strengths and preferentially in a co-planar configuration. However, in the case of catechol, the adsorption can also take place through a bidentate non-planar configuration that promptly dissociates, leading to a thermodynamically stable bidentate product.
The dissociation and transfer of the phenolic H ph to one of the framework O atoms connected to the Al is considered to be the first step of the tautomerization mechanism. The activation energy of this process ranges between 41 and 49 kJ/mol for phenol and between 12 and 55 kJ/mol for catechol, considering non-and co-planar adsorptions. The molecules tend to keep the co-planar configuration after the dissociation of the O ph H ph bond, which is followed by a rotation of the molecule over the Al centre. The rotation of the phenoxide intermediate requires between 0 and 12 kJ/mol, avoiding the regeneration of phenol and favouring the transfer of H ph to the C atoms in ortho position. The energy barrier of the H ph transfer to the C atom is only 22-27 kJ/mol for either phenol or catechol, considering both non-and co-planar adsorptions. In the case of catechol, the transfer is kinetically and thermodynamically more likely to involve the C atom that is not binding the second O 2ph . Thus, we have shown that the zeolite-mediated tau-tomerization of phenolic monomers is a viable process under mild temperatures, which changes the hybridization of the C atom from sp 2 to sp 3 , effectively disabling the aromatic conjugation. Finally, this work may help to understand better the role of the Lewis acid sites of zeolite when this material is used as a support during the hydroprocessing of lignin-derived compounds, and may contribute to better designed processes.