Surface reaction kinetics of the methanol synthesis and the water gas shift reaction on Cu/ ZnO/Al 2 O 3 †

A three-site mean-field extended microkinetic model was developed based on ab initio DFT calculations from the literature, in order to simulate the conversion of syngas (H 2 /CO/CO 2 ) to methanol on Cu (211) and Cu/Zn (211). The reaction network consists of 25 reversible reactions, including CO and CO 2 hydrogenation to methanol and the water-gas shift reaction. Catalyst structural changes are also considered in the model. Experiments were performed in a plug flow reactor on Cu/ZnO/Al 2 O 3 at various gas hourly space velocities (24 – 40 L h − 1 g cat − 1 ), temperatures (210 – 260 ° C), pressures (40 – 60 bar), hydrogen feed concentrations (35 – 60% v/v), CO feed concentrations (3 – 30% v/v), and CO 2 feed concentrations (0 – 20% v/v). These experiments, together with experimental data from the literature, were used for a broad validation of the model (a total of 690 points), which adequately reproduced the measurements. A degree of rate control analysis showed that the hydrogenation of formic acid is the major rate controlling step, and formate is the most sensitive surface species. The developed model contributes to the understanding of the reaction kinetics, and should be applicable for industrial processes ( e.g. scale-up and optimization).


Introduction
Methanol is a viable high-value energy carrier and a key intermediate in the chemical industry, as it can be used in many energy and chemical production applications. 1,2ethanol is synthesized from syngas (H 2 /CO/CO 2 ), which is industrially produced from natural gas and coal. 3However, because of concerns regarding CO 2 emissions and encouragements to the development of sustainable processes using renewable energy sources, alternatives to the generation of syngas are in the focus of worldwide research.Promising approaches include the gasification of biomass 4,5 and the use of green electricity (e.g. from solar and wind power) to generate hydrogen via water electrolysis, 6 further mixing it with carbon dioxide (e.g. from the flue gas of a process). 7hereby, the importance of the methanol synthesis will increase, as the conversion of hydrogen to secondary energy carriers is crucial. 8A scheme with the mentioned process steps is shown in Fig. 1.
2][13][14][15][16] Graaf et al. 11 proposed a model considering the three global reactions (eqn (1)-( 3)), one active site for CO and CO 2 adsorption and one for H 2 and H 2 O adsorption.Vanden Bussche and Froment 12 proposed a model considering CO 2 as the main carbon source of methanol and neglecting CO hydrogenation, but maintained the same active sites as Graaf et al. 11 Ovesen et al. 13 considered the three global reactions (eqn (1)-( 3)), one active site for CO adsorption and another one for CO 2 adsorption, and proposed structural changes in the catalyst active sites depending on the gas composition.Park et al. 14 considered the three global reactions (eqn (1)-( 3)) and three active sites (for CO adsorption, for CO 2 adsorption, and for H 2 /H 2 O adsorption).Seidel et al. 15 considered all three reactions and three active sites, applying the structural changes proposed by Ovesen et al. 13 as well.In a recent work, Slotboom et al. 16 proposed a model with a reduced amount of parameters (6)  to decrease identifiability problems.In their approach, CO hydrogenation (eqn (1)) was neglected, and three active sites were considered.
Although each model has its particular considerations, they all have rate determining steps (RDS) assumptions, which means that groups of parameters are lumped and fitted to experimental data.With this method, different effects may be merged with kinetics and cause the models to diverge outside the training region, which is often narrow.
Microkinetic modeling takes into account the chemistry behind the process, 17 being useful for better understanding of the system, for process optimization and for catalyst development.Based on data derived from first principles density functional theory (DFT) during the last decade, detailed microkinetic models have been proposed for methanol synthesis, [18][19][20][21][22][23] in which different surface reaction paths are considered and all reactions are potentially rate limiting.Grabow and Mavrikakis 18 proposed a mechanism for CO/CO 2 hydrogenation and the WGSR on Cu (111).In the adsorption steps, sticking coefficients equal to one were considered.The authors concluded that a more open surface (e.g.Cu (110), Cu (100), Cu (211)) could better represent the catalyst active area, and that the synergic effect of ZnO has to be taken into account, both conclusions being later demonstrated by Behrens et al. 24 Van Rensburg et al. 19 used previously reported DFT data 24,25 to test microkinetic mechanisms for CO and CO 2 hydrogenation on different surfaces (Cu (111), Cu (211), and Cu/Zn (211)), without including a mechanism for the WGSR.Liu et al. 20 compared the mechanisms on Cu 2 O (111) and on Cu (111), and concluded that CO hydrogenation is faster on Cu 2 O, whereas CO 2 hydrogenation is the dominating path on Cu (111).Park et al. 21proposed a mechanism on Cu (211), in which adsorption and activation energies were taken from DFTderived data, and pre-exponential factors were fitted to experimental data.Xu et al. 22 proposed a mechanism on Cu (211) and compared three situations: clean surface, preadsorbed O*, and preadsorbed OH*.The conclusion was that the preadsorbed species create faster reaction paths for the hydrogenation of formate.Huš et al. 23 developed a mechanism on Cu (111) for different metal oxides (Zn 3 O 3 , Cr 3 O 3 , Fe 3 O 3 , Mg 3 O 3 ) combined with copper, and concluded that the Zn 3 O 3 /Cu system has a superior performance.
7][28][29] To the best of our knowledge, only Huš et al. 23 validated their model with own experimental data.The other models cited above were either validated with experiments reported in literature, 11 or there was no validation at all.
Summarizing the modeling studies for methanol synthesis in the literature, to our knowledge, yet no published kinetic model exists that accounts for all these key features in methanol synthesis on Cu/ZnO/Al 2 O 3 : • The global mechanism considers reactions (1)-(3).
• All surface reactions are potentially rate limiting (detailed microkinetic approach).
• The facet of Cu (211), a more open and active surface of the catalyst, is considered.
• Three-site approach: besides the Cu (211) active site (site a), the synergy of zinc is taken into account with the assumption of a Cu/Zn or Cu/Zn δ+ (211) active center (site b), and a separate site is considered for the adsorption of H 2 and H 2 O (site c) to avoid the adsorption inhibition caused by formate.
• The structural changes of the catalyst, which are dependent on temperature and gas phase composition variations, are quantitatively taken into account.
• An extensive experimental validation is made covering the most important parameters in methanol synthesis: pressure, temperature, gas hourly space velocity (GHSV), and the feed composition (H 2 , CO and CO 2 ).
In this work, a multiscale kinetic model for the methanol synthesis and the WGSR is presented, in which all these aspects are considered.A wide experimental validation shall allow for simulating the methanol synthesis at different operating conditions.
The PFR set-up consisted of a stainless steel tube with 460 mm length, an inner diameter of 12 mm, and an inner concentric tube (2 mm) for temperature measurements in the axial direction.The feed gases were hydrogen (99.999% v/v), carbon monoxide (99.97% v/v), nitrogen (99.9999% v/v), and a mixture of carbon dioxide and nitrogen (50 : 50 ± 1.0% v/v) (Air Liquide Germany GmbH).The reactant gases supply was regulated via mass flow controllers (MFCs, Bronkhorst High Tech), by using proportional-integral-derivative (PID) control.The MFCs were calibrated with a flowmeter (Defender 530+, Mesalabs, standard error: 1.0% v/v).Both reactants (via bypass) and products were analyzed with a Fourier transform infrared spectrometer (FTIR, Gasmet CX4000).A flow diagram of the experimental setup is shown in Fig. S1.† The reactor was filled with a commercial CuO/ZnO/ Al 2 O 3 catalyst provided by an industrial partner, which was crushed and sieved to a particle size range between 250 and 500 μm.As the methanol synthesis is exothermic, and in order to avoid hot spots and to ensure isothermal operation, five portions of 0.30 g catalyst (in total 1.50 g), were separately mixed, each with 8.18 g of silicon carbide (SiC, Hausen Mineraliengroßhandel GmbH) (in total 40.90 g).Each mixture was then consecutively filled into the reactor, forming a catalytic bed length of 200 mm.Pure SiC completed the upper and lower ends of the bed.
The catalyst was activated as follows: a volume flow of 300 mL S min −1 containing 5% v/v of H 2 in N 2 was applied to the reactor, and the system was heated from 100 to 200 °C at a heating rate of 20 °C h −1 .This temperature was hold for one hour, followed by further heating to 240 °C at a heating rate of 12 °C h −1 .Finally, the H 2 concentration in the flow was increased to 50% v/v, maintaining the same total flow rate for one more hour.
In order to obtain a stable catalyst in steady-state conditions, the reactor was operated for 320 h at 40 bar, different temperatures (210-260 °C), and different feed gas compositions, before starting the measurements reported here.
The temperature axial profile was measured using a type-K thermocouple (NiCr-Ni) for the two data points with the highest methanol productivity, and thus the highest energy release due to the exothermic reactions.Since the maximum temperature difference was lower than 2 °C, the assumption of isothermal conditions is reasonable.The temperature profiles are presented in Fig. S2.† 3 Kinetic model development and numerical simulation

Three-site surface reaction mechanism
The developed kinetic model is based on the DFT calculations of Studt et al. 30,31 for CO and CO 2 hydrogenation to methanol, and the WGSR.The adsorption energies of CO 2 * and H 2 O* were taken from Polierer et al. 32 Studt et al. 30,31 applied the calculations on two stepped model surfaces: Cu (211), denoted in the presented mechanism as "site (a)", and a fully Zn-covered Cu (211), denoted "site (b)".The (211) facets were chosen, as it was found through experiments and theoretical calculations that they are the most active surfaces for the methanol synthesis. 18,22,24,30,31Since surface defects scale linearly with the overall observed activity, the reactivity of other facets (e.g.111) can be neglected. 24Besides, from a modeling point of view, the consideration of a most representative (i.e.most active) single facet is desirable, because an additional considered facet would double the amount of reactions and surface species, which makes the model more complex and requires a higher computational effort.
In this work, all reaction paths studied by Studt et al. 30,31 were originally implemented and tested for different operating conditions.The initial tests confirmed remarks proposed by Studt et al., 30 that CO hydrogenation on Cu/Zn (211) and CO 2 hydrogenation on Cu (211) are negligible.Therefore, CO hydrogenation on Cu/Zn (211) was eliminated from the final microkinetic model.The CO 2 hydrogenation on Cu (211), however, was kept in the model, as it describes the accumulation of formate on the Cu (211) surface, which reduces the number of free sites and could therefore slow down other reactions, e.g.CO hydrogenation.
Studt et al. 31 made DFT calculations for four WGS reaction pathways on both Cu (211) and Cu/Zn (211) surfaces: the redox mechanism, the water-assisted redox mechanism, the carboxyl mechanism, and the waterassisted carboxyl mechanism (totalizing eight possible reaction routes).After implementing and testing all these eight possible reaction pathways for different conditions, it was confirmed that the water-assisted carboxyl mechanism is dominant, and it is active on both Cu (211) and Cu/Zn (211).Thus, only this WGS reaction pathway (on both surfaces) is taken into account in the final kinetic model.After all, non-relevant reaction paths will only add complexity and increase the computational time of the simulations without contributing to the accuracy of the results.
It is known that formate (HCOO*) is able to cover a significant part of the catalyst surface, being an intermediate for CO 2 hydrogenation and inhibiting other reactions, such as CO hydrogenation. 30It is, however, unlikely that formate inhibits hydrogen adsorption, because there should be still enough small sites available Reaction Chemistry & Engineering Paper between adsorbed formate.A comparable situation is seen in the ammonia synthesis, in which it was shown that the interaction (or inhibition) of nitrogen with hydrogen is not significant. 33,34Another important feature is that, in the CO 2 hydrogenation, the decomposition of H 2 COOH* into H 2 CO* and OH* would need an additional free site, and could therefore be inhibited by high formate concentrations.On the other hand, this does not seem to be realistic, as H 2 COOH* is a large molecule that should not need extra space for its decomposition. 30Therefore, a third Cu (211) site is considered in our model, denoted as "site (c)", which is available for hydrogen and water adsorption.Similar approaches were published by other groups. 15,16,30he The reaction network of the carbon-containing species is shown in Fig. 2. It consists of five reaction pathways: CO hydrogenation on Cu (211), CO 2 hydrogenation on Cu (211), WGSR (water-assisted carboxyl mechanism) on Cu (211), CO 2 hydrogenation on Cu/Zn (211), and WGSR (water-assisted carboxyl mechanism) on Cu/Zn (211).

Kinetic equations
The catalyst surface is modeled considering a random distribution of the adsorbed species (mean-field approximation).A surface reaction is expressed as: With Here, N g and N s are the number of gaseous and surface species, respectively, χ i is the respective species i, v′ ik and v″ ik are the stoichiometric coefficients of the reactants and the Fig. 2 Reaction network of the carbon-containing species in the methanol synthesis and the WGSR.

Reaction Chemistry & Engineering Paper
products (species i in reaction k), respectively, and v ik is the stoichiometric gain of species i in reaction k.
As the methanol synthesis is typically operated at elevated pressures (50-100 bar), 9 the ideal gas consideration may give partial pressures that differ significantly from the actual fugacities of the gases.Slotboom et al. 16 reported deviations up to 10% comparing ideal and real gas approaches.Therefore, the Peng-Robinson equation of state 35 is used in our model to calculate the fugacities, using binary interaction parameters (k ij ) and other necessary data reported in literature, 36,37 and including an effective hydrogen acentric factor (ω = −0.05). 38n the present kinetic model, all surface species are considered to occupy a single site (σ i = 1 for all species).The surface coverage of a species i (θ i ) on a specific active site represents the fraction of this site that is occupied by that species.The calculation of the surface coverages is shown in eqn (6), and the sum of all coverages from a specific site must be 1 (eqn (7)).
where c i is the concentration of surface species i, σ i is the number of surface sites occupied by species i, and Γ is the surface site density.The turnover rate r (eqn (8)) consists in three multiplying functions: one dependent on temperature (F T ), 39,40 one dependent on the gaseous species fugacities (F G ), and one dependent on the surface species coverages (F S ).
Here, T is the reaction temperature, β is a correction due to the thermodynamic consistency (see section 3.3), k b is the Boltzmann constant, h is the Planck constant, E A is the reaction activation energy, ΔS ≠ is the reaction entropy barrier, f j is the fugacity of the gas component j, p 0 is the reference pressure (1 bar), R is the universal gas constant, and θ i is the surface coverage of species i. ϕ i represents the fraction of the site type of surface species i in relation to the total number of sites for carbon-containing compounds (sites a and b).Substituting eqn ( 9)-( 11) into eqn (12), the turnover rate of a reversible reaction k is: where the superscripts + and − refer to the forward and the reverse reaction, respectively.The reaction rate is related to the turnover rate by: The dependency of (ϕ i ) to the site type is shown in Table 1.The estimation of the zinc coverage (ϕ Zn ) is discussed in section 3.3.

Thermodynamic consistency
The microkinetic model has to correctly predict the thermodynamic equilibrium.The objective of the thermodynamic consistency corrections is to ensure reversibility of each elementary step according to the properties of the gas-phase species involved, which are known.The method described here is an adapted version of an approach developed in the Deutschmann's group, 41,42 considering a temperature operating range between T 1 = 200 °C and T 2 = 300 °C.
Goos et al. 43 reported thermodynamic data of gas species involved in the methanol synthesis.The free Gibbs energy function G 0 j of each gas species j at the reference pressure of 1 bar is: Here, a 1−7,j are compound-specific constants. 43The free Gibbs energy variation of the three global reactions (eqn ( 15)-( 17)) is: Table 1 Values of ϕ i depending on the site type

Site type
The assumption of constant heat capacity (c p ) was made for the global reactions in the temperature range of 200 to 300 °C.With this consideration, the Gibbs function was reduced from seven to three parameters (A 1−3 ) (eqn (18) and ( 19)) with the least square regression method (eqn (20)).
where The estimated regression parameters are summarized in Table S2.† When comparing the three-parameter functions with the seven-parameter ones, the average relative error of ΔG 0 m,3p was below 0.002%, and the maximum relative error was 0.007%.Therefore, a constant heat capacity sufficiently describes the free Gibbs energy change of the reactions involved in the methanol synthesis between 200 and 300 °C.
From DFT calculations, the free Gibbs energy change (ΔG 0 k, DFT ) of a reversible surface reaction k is given in the form: The free Gibbs energy change (ΔG 0 k,DFT ) of a global reaction pathway m (described in Fig. 2) is then calculated: Here, N r is the number of reactions, ζ k,m is the stoichiometric coefficient of a reversible surface reaction k in the global reaction pathway m.
In eqn ( 21) and ( 22) there is no term multiplying T•ln(T) for ΔG 0 m,DFT , like there is in eqn ( 18) and ( 19) for ΔG 0 m,3p .Each ΔG 0 m,DFT needs to be modified with the addition of parameter β multiplying T•ln(T), so that the equations are able to match.The calculation of β can also be seen in other thermodynamic consistency processes of surface kinetic mechanisms reported in literature. 41,42The new free Gibbs energy change (ΔG 0,TC m,DFT ) is calculated with the corrected terms ( ), as shown in eqn (23), and these terms are estimated so that eqn (24) holds for all reaction pathways: In eqn (24), the two functions will only be equal for a range of different temperatures (200-300 °C) if their corresponding terms match, namely the independent terms (eqn ( 25)), the terms accompanying T (eqn ( 26)), and terms accompanying T•ln(T) (eqn ( 27)).
where q 1−3,m represent the thermodynamic constraints.As in most microkinetic models, this is an underdetermined algebraic system, because there are 150 variables (6 parameters × 25 reactions) and only 15 equations (eqn ( 25)-( 27)) for the five reaction pathways).Herrera Delgado et al. 42 proposed an objective function that minimizes the individual corrections of E A , the pre-exponential factor (a term which contains ΔS ≠ ), and β.In this work, however, it was preferred to minimize the difference between the corrected Gibbs energy barrier (ΔG ≠,TC = E TC A − TΔS ≠TC − T•ln T•β TC ) and the original DFT-based one (G ≠,Orig = E Orig A − T•ΔS ≠,Orig ), for both forward and reverse reactions.This approach was chosen, because the model is more sensitive to modifications in the ΔG ≠ than in its individual parameters.Besides, the model is also sensitive to ΔG changes of all surface reactions (ΔG = ΔG ≠,+ − ΔG ≠,− ), even fast steps usually in equilibrium, as changing the ΔG will alter this equilibrium, which affects the whole mechanism.The constrained objective function is shown in eqn (28).
This journal is © The Royal Society of Chemistry 2021 Here, w k are selectable weights, which are chosen to protect the most sensitive reactions against changes.In this work, the weights of reactions R1, R2, R14, R16, and R17 were set to 10 5 and the other weights were set to one.This minimization problem can be solved with the method of the Lagrange multipliers (further explanation is given in the ESI †).
In Fig. 3, the capability of the model to predict the equilibrium is presented.The methanol concentration of the equilibrium for different operating conditions is calculated with Aspen Plus, using the RGibbs approach.Simulations with the same operating conditions are made with the microkinetic model considering a sufficiently long PFR, in order to achieve the equilibrium.When comparing the values, the conclusion is that the equilibrium is accurately predicted by the model, and the slight overestimations are probably due to rounding numbers and small differences in the thermodynamic data.

Estimation of the active sites distribution
A significant number of experimental observations have shown that the different active sites of the Cu/ZnO-based catalysts are adjusted dynamically to the operating conditions. 27,44It is therefore relevant to correctly model this phenomenon, in order to estimate the fraction of Cu (site a) and Cu/Zn (site b) on the surface for different operating conditions.After all, this quantification has a high impact on the simulation of the methanol synthesis.
The formation of a Cu-Zn alloy by reduction of zinc oxide and migration to the copper bulk can be described by the following reactions: 13,29 ZnO (s) + CO (g) ⇄ Zn (s) + CO 2(g) ΔG 0 (220 °C) = 60.89kJ mol −1 (29) Kuld et al. 29 proposed a detailed method to estimate the zinc fraction on the surface of the catalyst.First, the solubility of zinc in the Cu-bulk (X Zn ) is calculated considering the equilibrium of the zinc reduction via carbon monoxide (eqn ( 29)).The effect of the lower atom coordination in nanoparticles was considered for both, zinc oxide and copper.
where ΔG 0 Znred.(T) is the free Gibbs energy change in the zinc reduction via CO at the reference pressure (1 bar) and the reaction temperature T, γ Zn is the activity coefficient of zinc in Cu, a CO 2 and a CO are the activities of CO 2 and CO respectively,  ZnO and  Cu represent the respective surface energy of zinc oxide and copper, M ZnO and M Cu are the molar masses of zinc oxide and copper, d ZnO and d Cu are the crystallite diameter of zinc oxide and copper, ρ ZnO and ρ Cu are the density of zinc oxide and copper, and a ZnO is the activity of zinc oxide.The values of activity coefficient, surface energies and crystallite diameters were reported by Kuld et al. 29 The activity of the zinc oxide is assumed to be 1.Although not specifically noted by the authors, this assumption has probably been made considering that the typical quantity of metallic zinc in the copper bulk is not significant compared to the zinc oxide bulk.
In typical industrial methanol production, the WGSR is generally in equilibrium.Therefore, as the free Gibbs energy change of the reduction of zinc via hydrogen (eqn (30)) is significantly higher than via carbon monoxide (eqn ( 29)), the most probable way that hydrogen and water affect the reduction of zinc is by changing the a CO /a CO 2 ratio through the WGSR.Kuld et al. 29 proposed an effective a CO /a CO 2 ratio to be used in eqn (31) to account for the H 2 /H 2 O effect.
Here, K WGSR is the equilibrium constant of the WGSR.The segregation of metallic zinc from the Cu-bulk into the catalyst surface is then considered (eqn (33)).
Zn (Cu-bulk) + Cu (surf.)⇄ Zn (surf.)+ Cu (Cu-bulk) (33) Reaction Chemistry & Engineering Paper Kuld et al. 29 performed DFT calculations of this segregation on different facets, taking into account Zn-Zn interactions.For the facet Cu (211), the authors reported enthalpy variations (ΔH 0 seg ) of −27.01 kJ mol −1 for a Zn-free surface, −18.36 kJ mol −1 for a 0.333 Zn monolayer (ML), and −8.71 kJ mol −1 for a 0.667 Zn monolayer (ML).The ΔH 0 seg is then calculated (in kJ mol −1 ) as a function of the zinc coverage on the surface (ϕ Zn ): The entropy change (ΔS 0 seg ) of the segregation process on a Cu (211) facet was estimated to be 7.1 J mol −1 K −1 , and effects of Zn-Zn interactions in the entropy were neglected. 29The zinc coverage on the surface (ϕ Zn ) is calculated by solving eqn (35), in which it is considered that the segregation of zinc to the surface is in equilibrium.
The activities of the gases can be represented by their fugacities.The ratio ), represent the gas reducing power (GRP).Kuld et al. 29 showed that the described method agrees well with experimental data in a GRP range between 0 to 15, which roughly corresponds to a ϕ Zn range between 0 and 0.30 (for a Cu (211) facet), and a X Zn range between 0 and 0.003.However, the actual GRP of typical methanol synthesis operations may be higher than 15, especially at a CO-rich syngas condition, in which the amount of water is significantly low.For example, for a synthesis gas containing H 2 /CO/CO 2 /CH 3 OH/H 2 O/N 2 = 60/22/5/1/0.3/11.9%v/v, the GPR would be 880.In this case, it is assumed that the zinc concentration in the copper bulk may be more significant, causing a non-negligible reduction of the zinc oxide bulk, and thus affecting its activity.
Thus, the value for the ZnO activity coefficient is limited by the segregation of zinc to the copper bulk, and a linear relation between a ZnO and X Zn is made (eqn (36)), in order to extend the application range of the method from Kuld et al. 29 for higher GRP.
Here, α express this linear relation, and a value of α = 1.5 is chosen, which roughly correspond to a maximum zinc coverage of 0.90, a limit value proposed elsewhere. 15ubstituting eqn (36) into eqn (31), representing the gas activities by their fugacities, and reorganizing the equation, the fraction of zinc in the copper bulk (X Zn ) is: In Fig. 4, the original and the modified functions are compared for different GRP at the temperatures of 220 and 250 °C.In Fig. 4A and C, the solubility of zinc in the Cu-bulk is shown, while in Fig. 4B and D the zinc coverage is presented.It is show that the modified function attenuate X Zn and ϕ Zn at a high GRP, while not influencing X Zn and ϕ Zn significantly at low GRP.
Finally, at CO 2 -rich operating conditions, the reverse water-gas shift reaction (rWGSR) is normally faster than the WGSR, and is usually not in equilibrium.In this case, the contribution of the direct reduction of ZnO by H 2 is probably higher, and a more complex analysis must be made. 29In order to keep using the method proposed by Kuld et al., 29 it was chosen to set a minimum zinc coverage of 0.30 to avoid underestimations of zinc coverage by these conditions.

Estimation of the active catalytic area
Because of the dynamic behavior of the Cu/ZnO/Al 2 O 3 and the three-site approach, characterization tests (e.g.N 2 O chemisorption) can just estimate an initial active catalytic area, since it changes depending on the experimental conditions.Still, these estimations serve as a reference in comparing different catalysts.The determination of the surface site density is a challenge for the same reason.Therefore, it was chosen to use the experimental data to estimate a specific catalyst site quantity (n M,Cat ) in terms of mol of active sites per catalyst mass unit (eqn (38)).The total active surface area (A Cat ) is then calculated by eqn (39).
Here, m Cat is the mass of the catalyst inside the reactor.By estimating the n M,Cat , the need to quantify the surface site density (Γ) is avoided, as shown in the next section.
The n M,Cat is estimated by minimizing the prediction errors of the methanol output with the experimental data.When the experimental data has values differ significantly (e.g.0.08% v/v and 12.00% v/v), a better distribution of each point's importance can be made with introducing weights.Common approaches are the inverse of squared experimental value, 45,46 the inverse of the squared simulated value, 47 or the inverse of experimental multiplied by simulated value. 48ere, the inverse of the squared experimental value was used as weights.The function fminsearch from Matlab was used, and the objective function is shown as follows. (37)

Reaction Chemistry & Engineering Paper
This journal is © The Royal Society of Chemistry 2021 Here, y n CH3OH;out and ŷn CH3OH;out are the experimental and simulated value of point n, respectively.

Reactor equations
In this work, the microkinetic model is applied to simulate steady-state operation of two types of reactor: a fixed-bed tube reactor (own experiments and literature data 14,16 ) and a CSTR (literature data 15 ).Isothermal operation was considered in both cases.
3.6.1 Fixed-bed tube reactor.In the tube reactor model, only variations along the reactor length are assumed, given the ratio between the diameter of the reactor and the particle size (24 ≤ d R /d p ≤ 48).The influence of back-mixing is neglected (plug flow reactor assumption, PFR).A total molar balance along the catalyst bed length (L) is calculated: where n ˙is the total gas mole flow, z is the axial direction, N r is the number of reactions.Substituting eqn ( 13) and ( 39) in eqn ( 41): The axial reactor profile of the molar fraction of each gaseous species j (y j ) is calculated via a component balance of the gas phase.
The coverage θ i of each surface species i at a certain point in time is calculated via a component balance of the surface.
Comments on solving this system of differential equations are given in the ESI.† 3.6.2Continuous stirred tank reactor.In the CSTR model, a total molar balance in the reactor is calculated: where dn/dt is the total mole accumulation in time, n ˙in is the mole flow entering the reactor, and n ˙out is the mole flow leaving the reactor.Assuming no gas accumulation in the reactor: The component mole balance in the reactor is calculated: where dy j /dt is the change in time of the mole fraction of component j, n is the total mole quantity in the gas phase, y j,in is the mole fraction of component j entering the reactor, and y j is the mole fraction of component j in the reactor.The mole quantity in the gas phase can be calculated with the Peng-Robinson equation of state. 35ike in the PFR, the coverages of the surface species are calculated by: Comments on solving this system of differential equations are given in the ESI.†

Sensitivity analysis
In order to evaluate the most sensitive reaction rate parameters in the kinetic model, the Campbell degree of rate control (DRC) method was applied. 49,50This method consists in slightly changing the Gibbs energy (G i ) of a surface intermediate or a transition state, while keeping the Gibbs energy of the other species G w≠i constant, and it has the advantage of maintaining the thermodynamic consistency of the model.For a set of reversible reactions, the degree of rate control of surface species or a transition state i (DRC i ) is defined as: Here, (r 5 + r 6 ) is the methanol production rate, and G 0 i is the free Gibbs energy of species i at the reference pressure (1 bar).The method of finite differences is used as an approximation to solve eqn (50), and a step δ = 0.01 kJ mol −1 was chosen.Eqn ( 51) is used to calculate the sensitivity of methanol generation, and eqn ( 52) is used for sensitivity of CO generation, which makes sense at high CO 2 content.
4 Results and discussion

Model validation
A microkinetic model for syngas (H 2 /CO/CO 2 ) conversion to methanol including a three-site approach and structural changes was successfully developed.The complete set of 25 reversible reactions and their respective parameters to calculate the turnover rates are summarized in Table 2.The thermodynamic consistency of the model was ensured (see section 3.3).The estimated value of the catalytic site quantity (n M,cat ) is 2.00 mol kg cat.−1 .
5][16] The operating conditions in the respective setups are significantly different from ours, which contributes to a broader validation range.In Table 3, the operating conditions of each setup is summarized.For the data from Seidel et al. 15 with H 2 /CO in feed, our model was initially overestimating methanol formation.Feeds with only H 2 /CO are an extreme case, when looking at the zinc coverage.As the catalyst used is not exactly the same, we assume that in the experiments of Seidel et al. 15 zinc is covering more surface, leaving less Cu sites for the CO hydrogenation.By adjusting the value of α to 1.17 (see eqn (36)) in this specific case, roughly corresponding to a maximum of 95% zinc coverage, the measurements from Seidel et al. 15 with H 2 /CO were reasonably predicted.
In Fig. 5, the normalized residues of the original model are shown for each carbon-containing species.By comparing simulations with experimental data, it was found that the WGSR is adequately predicted.However, if the operating conditions favor rWGSR, generally for feed ratios of CO 2 /CO X (y ¯CO 2 ,0 ) higher than 0.65 (CO X = CO + CO 2 ), the simulation of CO production through the rWGSR was significantly higher than the experimental values, as shown in Fig. 5A.This suggests that the Gibbs energy barrier (ΔG ≠ = E A − T•ΔS ≠ ) is influenced by higher concentrations of CO 2 and H 2 O or by surface intermediates derived from them, namely HCOO* and OH*, respectively.Different approaches were tested to improve the simulations in this region, including the addition of coverage dependency terms.The solution found was to add 15.44 kJ mol −1 to the activated complex energies of the most sensitive reactions of the rWGSR (R24 and R25), if the feed ratio of CO 2 /CO X is higher than 0.65 (see Table 2).With this procedure, the model remains thermodynamically consistent.
In Fig. 6, the normalized residues of the two-case model are shown for each carbon-containing species.The model quantitatively reproduces the 690 measurements, with the majority of the data points being inside the ±20% lines (58% of the methanol points and ≥90% of CO and CO 2 points).CO 2 simulation error is always smaller than 40%, and only 3.4% of CO data points and 17% of methanol data points are outside the ± 40% lines.The sum of the relative squared errors (χ 2 ), the mean error (ME j ) and the mean squared error of the predictions are calculated for CO (j = 1), CO 2 (j = 2), and methanol (j = 3), and are shown in Table 4.
Considering all experiments, the mean error (ME) values of CO, CO 2 and methanol concentration are 8.0, 3.3, and 21.7%, respectively.These values are even lower when looking only to mixed feeds (H 2 /CO/CO 2 ).The mean squared errors (MSE) are also significantly low (see Table 4).
In Fig. 7, it is shown the experimental and simulated values of methanol output concentration for different conditions and setups.It can be seen that the simulations are significantly close to the experiments and the trends are

Forward reaction
Reverse reaction Reaction Chemistry & Engineering Paper adequately predicted.In Fig. 8, the error of the prediction of the carbon-containing compounds is shown as a function of y ¯CO 2 ,0 (x-axis) and temperature (y-axis) for our experiments at the operating conditions of 40 bar, 24 L h −1 g cat −1 and H 2 /CO x /N 2 ≈ 45.3/14.3/40.4% v/v.The model simulates CO and CO 2 accurately for the entire studied region.The simulation of methanol is also reasonable, with slight underestimations at high temperatures and low y ¯CO 2 ,0 , and moderate underestimations at low temperatures and high y ¯CO 2 ,0 .This leads to the conclusion that there is some positive effect on CO 2 hydrogenation at low temperature when the concentration of CO 2 is increased, which is not reflected in this model.This effect is shown in both our own experiments and experiments from literature. 51Still, the model adequately simulates this region (low temperature and high y ¯CO 2 ,0 ) for lower values of GHSV is lower (which means higher conversion).This is demonstrated in the simulation of the experiments from Slotboom et al. 16 which have no CO in feed (y ¯CO 2 ,0 = 1) and are reasonably reproduced.

Reaction flow and sensitivity analysis
The validated model was used to simulate the methanol synthesis at an extended range of conditions.In Fig. 9, it is shown the turnover frequency of the different reaction paths and the CO X conversion along the reactor.The operating conditions are 60 bar, 220 °C, 4.8 L S h −1 g cat −1 , and a feed concentration of H 2 /CO x = 80/20% v/v.Fig. 9 is complemented by the coverages of the surface species (Fig. 10) and the zinc coverage along the reactor (Fig. 11).In Fig. S6-S8 † analogous diagrams are shown for an operating temperature of 250 °C.
According to the simulations, CO hydrogenation on site a (Cu) is only relevant at CO-rich feeds (e.g.20% of the methanol production at y ¯CO 2 ,0 = 0.25), CO 2 hydrogenation on site a (Cu) does not occur significantly at any condition, and CO 2 hydrogenation on site b (Cu/Zn) is the main reaction path for the production of methanol, in agreement with findings from DFT studies. 24,30t CO-rich conditions (e.g.y ¯CO 2 ,0 ≤ 0.50), CO 2 conversion and (consequently) water generation are fast within the initial 10 cm of the reactor, mainly due to CO 2 hydrogenation, but probably also due to contributing rWGSR (depending on the CO 2 /CO x ratio and the temperature).Furthermore, the WGSR rate increases rapidly along the reactor, and after a certain amount of water has been produced (it is dependent on the operating conditions), the WGSR rate is the approximately equal to the CO 2 hydrogenation rate.From this axial position on (ca. 10 cm, Fig. 9A and B), CO 2 concentration remains constant, because it is consumed in the CO 2 hydrogenation but regenerated in the WGSR.Therefore, from this axial position on, only CO is converted, through the combination of WGSR and CO 2 hydrogenation and through direct CO hydrogenation.
With more CO 2 content in the feed (e.g.y ¯CO 2 ,0 ≥ 0.75), a higher conversion of CO 2 is achieved, which implies that higher amounts of water are generated, and the methanol synthesis is severely slowed down.At CO 2 -rich feeds, there is significant CO production via the rWGSR, increasing with temperature (see Fig. S7 †).
At constant temperature and pressure, the decrease of CO 2 hydrogenation rate on Cu/Zn, along the reactor length, has three main causes: • The decrease of H (c) coverage, which is caused by H 2(g) consumption; • The zinc coverage (ϕ Zn ) decrease, caused by the generation of water and consumption of H 2 and CO, 27,29 which significantly changes the gas reducing power (eqn (32)).
• The product inhibition with the increase of H 3 CO (b) and OH (b) , because of methanol and water accumulation, respectively.
From our simulations, the product inhibition has the highest effect in reducing the reaction rate.This is particularly relevant for CO 2 -rich feeds, in which a significant water accumulation takes place and, thus, OH (c) reduces the amount of free sites (c), as shown in Fig. 10C and D.
In a recent work, 52 a methanol-assisted autocatalytic mechanism was proposed, with water as the only effective inhibitor of the methanol synthesis.If appropriate DFT calculations were accessible, our model could be extended accordingly and this possibility could be investigated further.
In Fig. 12, the positive effects of pressure and CO 2 concentration, and the negative effect of temperature on formate coverage on Cu and on Cu/Zn are shown.Formate requires lower temperature, higher pressure, and higher CO 2 concentration to block the majority of Cu sites (Fig. 12A and B), but it covers most of the Cu/Zn sites even at mild conditions (Fig. 12C and D).Therefore, moderate consumption of CO 2 changes the formate coverage on Cu/Zn only slightly, and therefore the reaction rate of CO 2 hydrogenation is not as much affected by CO 2 consumption as by the factors previously discussed.This was experimentally demonstrated by varying the feed concentration of CO 2 while maintaining the other conditions constant, as shown in Fig. 13.

Reaction Chemistry & Engineering Paper
The method of degree of rate control (DRC) is applied to investigate the sensitivity of the methanol production in relation to the free Gibbs energy of each surface intermediate and of each transition state of the reversible reactions.This analysis is shown at 60 bar (Fig. 14) and different temperatures considering a gas phase concentration at low conversion: H 2 /CO x /CH 3 OH/H 2 O = 79.8/19.8/0.2/0.2%v/v.A separate analysis is made for the case of little CO content in the gas (Fig. 15), which includes the sensitivity of both CO and methanol generation.The DRC was also applied to a   2) for the entire operating region under study, which is in agreement with other DFT-derived models 19,53 and lumped kinetic models. 11,16Xu et al. 22 concluded that the hydrogenation of formic acid is the RDS for the CO 2 hydrogenation at lower formate coverages, and the H 2 COOH* association (R19, Table 2) is the RDS for high formate coverages, the latter conclusion being a result of the one-site approach made by the authors (see section 3.1).Grabow et al. 18 and Park et al. 21found the hydrogenation of H 2 CO* to be the RDS on Cu (111) without considering the zinc influence.Finally, the model presented here shows that the hydrogenation of formate (R15), a typically assumed RDS in formal kinetic models, 12,15 has a reasonable sensitivity in our model (0.20-0.35), but still far behind the sensitivity of R17 (0.50-0.70).
In the CO hydrogenation on Cu, the most sensible reaction of this microkinetic model is the hydrogenation of HCO (a) (R11), which is in agreement with findings of Van Rensburg et al. 19 In formal kinetic modeling studies, Graaf et al. 11 proposed the hydrogenation of H 2 CO* (R20) to be the RDS, while Seidel et al. 15 assumed that the RDS is the hydrogenation of H 3 CO* (R6).
From our DRC analysis of CO generation sensitivity at CO 2 /CO x = 0.987 (Fig. 15), the reaction between CO 2 * and H 2 O* on Cu (R24) is the RDS of the rWGSR.Other reported models that include the WGSR usually consider the redox mechanism 12,15 or the carboxyl mechanism without water assistance 16 to be the RDS, which have been tested here with DFT data and also evaluated elsewhere as not significant if compared to the water-assisted carboxyl mechanism. 31egarding the participating intermediates, our findings suggest that formate on Cu/Zn (HCOO (b) ) is the most sensitive species, as it is the most abundant species, and participates on the second most sensitive reaction (formate hydrogenation), which produces the formic acid (HCOOH (b) ).
The rWGSR is mostly sensitive to formate adsorbed on Cu, although formate itself is a spectator of this reaction path, suggesting that formate may inhibit the rWGSR at a certain level.This hypothesis is further supported by the argument

Summary and conclusions
A thermodynamically consistent microkinetic model was successfully developed, based on first principles DFT derived kinetic data for the methanol synthesis and the WGSR on Cu (211) and Cu/Zn (211).A method for the quantification of the zinc coverage dependent on the operating conditions, proposed in the literature, is applied here, and is slightly modified to improve the estimation for gas compositions rich in H 2 and CO.
With an extensive validation consisting in 359 own experiments (available in the ESI †) and experiments from three different sources in literature (i.e.139 data points from Seidel et al., 15 98 data points from Park et al., 14 and 94 data points from Slotboom et al. 16 ), the model reproduces the system quantitatively in a broad range of relevant conditions, showing discrepancies only for the combination of low temperature, high CO 2 /CO x concentration in feed and high GHSV.The proposed model is based on theoretical calculations, and we believe it has a high change of accurately predicting the methanol synthesis outside the validation region.
The reaction flow analysis showed that methanol is mainly formed from CO 2 hydrogenation on site b (Cu/Zn), and that CO conversion is mostly due to the WGSR on site a (Cu).At CO-rich conditions, direct CO hydrogenation is responsible   for some of the methanol generation (e.g.20% of the methanol production for CO 2 /CO X = 0.25).On the other hand, this reaction pathway is strongly inhibited by formate accumulation on Cu surface at higher CO 2 concentrations.The model suggests that formation of methanol and water leads to an accumulation of H 3 CO (b) and OH (c) , respectively, with both slowing down the overall reaction.This is particularly significant in the case of CO 2 -rich feeds, as high amounts of water are generated in the process.We assume that the productivity should be increased by using reactor designs that enable product extraction in situ, such as using membranes, or integrating reaction and product condensation steps.This is especially encouraged for the conversion of CO 2 to methanol.
The sensitivity analysis, using the method of the degree of rate control (DRC), pointed out that the formic acid hydrogenation on site b (Cu/Zn) (HCOOH (b) ) is the most sensitive step (DRC around 0.60).Formate (HCOO (a) ) was found to be the most sensitive intermediate in the rWGSR, although it does not participate in this reaction, which suggests that the rWGSR is slowed down due to formate blocking of free sites.With this finding and based on the premise that formate only reacts further to produce methanol, we formulate the following hypothesis: If a modification on the Cu/Zn-based catalysts is realized, in which formate binds stronger on the Cu site and achieves coverages closer to 1, the rWGSR might be more effectively inhibited, and CO 2 conversion at CO 2 -rich systems should be enhanced.
The presented microkinetic model can be further extended with additional reactions paths if DFT calculations are available (e.g.methanol dehydration to DME or methanol-assisted autocatalytic reaction paths).
The quantification of the zinc coverage may be improved by an approach on a higher level of theory for the zinc oxide activity, by an operando analysis to check which is the maximum zinc coverage on the surface at a H 2 /CO feed, and with an approach for non-equilibrated WGSR and domination of zinc reduction via hydrogen (e.g. for high CO 2 / CO X ratio in feed).
ΔG 0 m,3p is the Gibbs energy change of the global reaction m considering three parameters, ΔH 0 Tr,m and ΔS 0 Tr, m are the enthalpy and entropy change of the global reaction m at the standard temperature (T r = 298.15K), Δc p,m is the heat capacity change of the global reaction m, and A 1−3,m are the regression parameters 1-3 of the global reaction m.

Fig. 3
Fig. 3 Methanol equilibrium concentration calculated with Aspen Plus and with the microkinetic model.Operating conditions: 60 bar, and a feed concentration of H 2 /CO x = 80/20% v/v.

Fig. 4
Fig. 4 Solubility of zinc in the Cu-bulk (A and C) and zinc coverage (B and D) as functions of the gas reducing power.

Fig. 7
Fig.7Experimental and simulated values of methanol output concentration at different conditions.Databases: A) this work.B) Seidel et al.15

Fig. 8
Fig. 8 Simulation error (color) as a function of CO 2 /CO x in feed and temperature.Operating conditions: 41 bar, GHSV = 24 L S h −1 g cat −1 ,

Fig. 11
Fig. 11 Zinc coverage along the methanol synthesis reactor with a length of 100 cm.Operating conditions: 220 °C, 60 bar, GHSV = 4.8 L S h −1 g cat

k
Coefficient n of the Gibbs term of reaction path m (eqn (19)) a n,jCoefficient n of the Gibbs term of gas species j (eqn (14)) a j Activity of gas component j Δc p,m Heat capacity change of global reaction m (kJ mol −1 K −1 ) d i Crystallite diameter of species i (m) DRC i Degree of rate control of a transition state or surface intermediate i E A,k Activation energy of reaction k (kJ mol −1 ) f j Fugacity of gas component j (Pa) G 0 j (T) Free Gibbs energy at 1 bar and temperature T (kJ mol −1 ) ΔG 0 R (T) Free Gibbs energy of (global/surface) reaction R at 1 bar and temperature T (kJ mol −1 ) ΔH 0 298.15KStandard reaction enthalpy (kJ mol −1 ) ΔH 0 seg Segregation enthalpy at 1 bar (kJ mol −1 ) (see eqn (33) and (34)) h Planck constant (6.62607 × 10 −34 J s) K R Equilibrium constant of reversible reaction R k b Boltzmann constant (1.38065 × 10 −23 J K −1 ) L Catalyst bed length (m) M i Molar mass of species i (kg mol −1 ) m Cat Catalyst mass (kg) N g Number of gas components (-) N r Number of surface reactions (-) N s Number of surface intermediates (-) n M,Cat Specific catalyst site quantity (mol kg −1 ) n ˙Total mole flow (mol s −1 ) p 0 Reference pressure (1 × 10 5 Pa) q n,m Thermodynamic constraint n of reaction pathway m R Universal gas constant (8.31446 × 10 −3 kJ mol −1 K −1 ) r k Turnover rate of reversible reaction k (s −1 ) ΔS ≠ Entropy barrier of reaction k (kJ mol −1 K −1 ) s ˙k Rate of reaction k (mol m −2 s −1 ) T Temperature (K) w k Selectable weight of reaction k (-) y j Mole fraction of gas component j (mol/mol) y n j,out Experimental value (point n) of mole fraction of gas component j of point n in the reactor outlet (mol/mol) y ˆn j,out Simulated value (point n) of mole fraction of gas component j of point n in the reactor outlet (mol/mol) X Zn Solubility of zinc in the copper bulk (mol/mol) Greek α Proposed linear relation parameter between the zinc oxide activity and the zinc solubility in the Cu bulk (-) β k Correction term of reaction k because of the thermodynamic consistency (-) Γ Surface site density (mol m −2 ) δ Step size of the finite differences method (kJ mol −1 ) γ Zn Activity coefficient of metallic zinc (-)  i Surface energy of species i (kJ m −2 ) ξ k,m Stoichiometric coefficient of a reversible surface reaction k in the global reaction pathway m θ i Surface coverage of species i (-) ρ i Density of species i (kg m 3 ) σ i Number of surface sites occupied by species i (-) ν′ ik Stoichiometric coefficient of reactant i in reaction k (-) ν″ ik Stoichiometric coefficient of product i in reaction k (-)

Table 4
Statistical indicators of the model performance in predicting the carbon-containing compounds Reaction Chemistry & EngineeringPaperν ik Stoichiometric gain of species i in reaction k (-) ϕ iFraction of the site type of surface species i in relation to the total number of sites for carboncontaining compounds (sites a and b) χ i Species i