Understanding catalytic CO 2 and CO conversion into methanol using computational ﬂ uid dynamics †

The kinetics of methanol synthesis from a mixture of CO 2 /CO/H 2 have been widely studied in the literature. Yet the role of direct CO hydrogenation is still unclear, in terms of predicting and developing an accurate kinetic model. To investigate, a computational ﬂ uid dynamics model has been developed, incorporating two distinct kinetic models, one which includes CO hydrogenation and one which does not. Including CO hydrogenation in the kinetic model provides a more complex interaction between the three involved reactions and can better predict potential inhibitions caused by the presence of H 2 O. This, however, increases the complexity of the kinetic model. The bene ﬁ t of applying a ﬂ uid dynamics model to study ﬁ xed bed reactors is demonstrated, as it o ﬀ ers unique insights into the spatial species concentration, temperature variations, and reaction rate magnitudes. The validated model is shown to be a powerful interrogative tool, capable of supporting system optimization across the catalyst and reactor engineering sectors.


Introduction
Increasing greenhouse gas emissions continually motivate research into offsetting anthropogenic emissions. Carbon Capture and Utilisation (CCU) has great potential in the long-term mitigation of CO 2 emissions, 1 especially when applied to the large scale production of fuels or bulk chemicals, such as the methanol (MeOH) economy concept, outlined by Olah et al. 2 Here, captured CO 2 can act as a MeOH feedstock, either through its direct hydrogenation, eqn (1), or through a two-step process where CO is initially produced through the Reverse Water-Gas Shi (RWGS) reaction, eqn (2), followed by subsequent CO hydrogenation to MeOH, eqn (3). By nature, the forward rate of all three reactions is exothermic, while their reverse rate is endothermic. R 1 : CO 2 + 3H 2 4 CH 3 OH + H 2 O, À49.2 kJ mol À1 (1) R 2 : CO + H 2 O 4 CO 2 + H 2 , À41.2 kJ mol À1 (2) R 3 : CO + 2H 2 4 CH 3 OH, À90.8 kJ mol À1 (3) impactful than CO 2 hydrogenation, can still occur. Its signicance, however, varies greatly in the literature; as an example, Grabow and Mavrikakis 34 predicted $1/3 of MeOH being produced from CO, while Sahibzada et al. 23 predicted that MeOH synthesis from CO 2 may be 20 times higher than MeOH synthesis from CO. Grabow and Mavrikakis 34 investigated the Cu(111) surface of the commercial catalyst, while in the DFT study of Higham et al., 37 the Cu(110) and the Cu(100) facets were investigated instead; it was identied that the Cu(110) and Cu(100) facets could be even more active for CO 2 dissociation into CO and subsequent hydrogenation to MeOH 37 compared to the Cu(111) facet, yet the impact of CO hydrogenation in MeOH synthesis was not dened. Wilkinson et al. 43 investigated the species concentration at different lengths along the catalytic bed, by having different catalyst loads in the reactor each time. The conclusions seem to indicate that the Forward Water-Gas Shi (FWGS) reaction followed by MeOH synthesis is a more favourable pathway for CO compared to direct hydrogenation. It is easy to see from these reports that to date, there is still no denitive consensus on how impactful CO hydrogenation is in MeOH synthesis. Another uncertainty concerned the active catalytic site for CO and CO 2 adsorption. Graaf et al. 7,44 introduced a KM based on a dual-site Langmuir-Hinshelwood mechanism, where CO and CO 2 adsorb competitively on the rst site and H 2 and H 2 O adsorb competitively on the second site. Graaf's KM has since become one of the most widely used models for MeOH synthesis kinetic investigations, 6,12,13,25,38,45 even though it is known to underestimate experimental MeOH synthesis rates. 6 However, different studies proposed that CO and CO 2 follow a non-competitive route to MeOH, with CO adsorbing on Cu 1+ sites and CO 2 adsorbing on Cu 0 sites. 8,23 Based on this, Park et al. 46 proposed an alternative KM, as shown in eqn (4)-(6) which refer to reactions R 1 -R 3 in eqn (1)-(3), respectively. In this KM, CO 2 and CO hydrogenation rates (eqn (4) and (6), respectively) are independent from one another. To dene the kinetic constants, the authors ran an extensive experimental session with the commercial Cu/ZnO/ Al 2 O 3 catalyst in a xed bed reactor, by varying the CO/CO 2 /H 2 feed composition, total ow rate, operating temperature and pressure.
Kinetic model from Park et al.: 46 Park's KM was validated in a subsequent paper, where it was integrated into a CFD model, with the focus on investigating the thermal behaviour of a xed bed reactor. 47 As a result of these ndings, the authors proposed an alternate bayonettype reactor be used, which demonstrated an enhanced thermal prole, reducing the peak temperature reached in the bed, and achieving a higher MeOH conversion. Their study demonstrated how CFD can guide experimental procedure. In another recent paper, a hybrid Cu/ZnO/Al 2 O 3 -ferrierite catalyst, used for direct DME synthesis from syngas via a MeOH intermediate, was investigated. 48 The study incorporated Park's KM into a Matlab code which investigated the mechanistic implications of the hybrid catalyst, thus identifying MeOH synthesis as the primary rate controlling step whilst complementing the experimental study. This required the kinetic constants of Park's reaction rates to be reevaluated for this new hybrid catalyst.
Another widely studied KM, [49][50][51] proposed by Bussche and Froment, 52 assumed all reactions to occur only on the Cu catalytic sites. Here, CO hydrogenation is omitted, as MeOH is assumed to be exclusively synthesised from CO 2 . In this KM, CO participates only through the FWGS reaction, serving as a CO 2 precursor for subsequent conversion to MeOH.
To identify which of the available models is better suited to predicting experimental results, several studies have compared them. Nestler et al. 53 and Wilkinson et al. 43 both compared the Bussche and Graaf KMs, but did not consider Park's KM. 46 Both concluded that Bussche's KM is more accurate, as the kinetic constants of Graaf's model were derived for an older, less active version of the commercial MeOH synthesis catalyst. 53 Following this, Nestler et al. 53 insisted that these KMs are only suitable for the specic range of pressures, stoichiometric numbers and CO/CO 2 ratios from which they were derived. They should not be extrapolated outside this range to minimise the simulation errors. Here, an alternative KM was also proposed using eqn (7) and (8), by applying a parameter tting to the experimental data set produced by Park et al. 46 Nestler's KM was established to be applicable to a wide range of industrially relevant MeOH synthesis operating conditions, specically a pressure range of 50-80 bar, temperature range of 473-593 K, and a varying range of CO 2 , CO, and H 2 feed compositions, including pure CO 2 feeds with no initial CO content. 53 On rst sight, the main difference between Park's KM (eqn (4)-(6)) and Nestler's KM (eqn (7) and (8)) is the inclusion of CO hydrogenation. Yet, they also have some additional mechanistic differences. The denominator of both KMs represents the adsorption of species on the active sites; as previously stated, Park's KM, eqn (4)-(6), considers CO 2 and CO to have two distinct sites, as the adsorption of one is independent from the adsorption of the other, i.e. the denominator of r CO 2 depends only on CO 2 while the denominator of r CO depends only on CO. Nestler's KM, eqn (7) and (8), considers a common catalytic adsorption site for CO and CO 2 , as the denominator of r CO 2 depends on the fugacity, here assumed to be equal to the partial pressure, of both CO 2 and CO. Both models consider a common site for H 2 and H 2 O. The reaction rate of the RWGS reaction is dependent only on the adsorption of CO 2 in Park's KM, while Nestler's KM considers both CO and CO 2 . Overall, Park's KM considers three distinct active sites, while Nestler's KM considers two.
As previously mentioned, CFD models can be used as an investigative tool, since their predicted ow prole can indicate the trends and the mechanisms of chemical reactions. In our previous work, one such CFD model was applied to study ethanol dehydration to ethylene on a xed bed reactor using a SAPO-34 catalyst. 54 By combining experimental and CFD studies, it was possible to identify the preferable mechanistic pathway of ethanol dehydration. Here, the previously presented CFD model will be adapted and applied to investigate the kinetic mechanisms of MeOH synthesis, using the KMs of Park and Nestler. To compare and validate the two KMs, in terms of species interactions and shis in the chemical equilibrium, the experimental results previously presented by Park et al. 46 will be used. The validation in this paper will focus mainly on the mechanistic implications of changing the composition of the feed gas, rather than the operating temperature, ow rate, and pressure. Thus, only the experimental runs with temperature of 523 K, ow rate of 8000 mL g cat À1 h À1 , and pressure of 50 bar are considered. Noticeably, both models were derived by tting their parameters to the same experimental setup and results, those of Park et al., 46 yet they assume different reaction mechanisms and roles for CO. Identication of the strengths and internal inconsistencies in these mechanisms can provide key information regarding the impact of CO hydrogenation in MeOH synthesis.

Computational model
In this section, the computational geometry, which is based on the experimental setup of Park et al., 46 will be presented. In addition, the methodology used to implement diffusion transfer and the KMs of Park and Nestler are analysed. Details of the computational theory and the setup of the porous medium, the boundary conditions, and the solution setup are presented in the ESI. †

Computational geometry
The Ansys® Academic Research, release 2019 R2, simulation soware was used for this study, with Fluent as the computational code. A 3D cylindrical geometry was designed with the same dimensions as the xed bed used in the experimental setup of Park et al., 46 where 0.4 g of Cu/ZnO/Al 2 O 3 catalyst (assumed density 1775 kg m À3 ) 52 was loaded in a 7 mm inner-diameter reactor tube. The estimated height for the cylindrical packed bed, based on the aforementioned density and the estimated volumetric porosity (see ESI †), is 15.6 mm. In the experimental setup, the catalyst was diluted with 1.2 g of inert Al 2 O 3 as a means to limit the exothermicity of the reaction. These inert species will also affect the ow dynamics and overall pressure drop along the packed bed. However, implementing the inert Al 2 O 3 diluent into the current porous medium model would increase the complexity of the code, since inert species should be dened as unreactive computational elements, homogeneously mixed with reactive elements. Thus, in the considered geometry, these inert species were not considered. The generated mesh consists of cubic elements of 0.1 mm in size, resulting in a total of 627k elements. The schematics of the reactor, along with the generated mesh, are presented in Fig. A1 of the ESI. † In addition, a mesh independency study is also presented in the ESI. † This study, apart from justifying the decision for the chosen mesh element size, also indicates the renement necessary to accurately capture the ow characteristics while minimizing errors and/or computational time. A steady-state simulation was considered, with laminar conditions assumed, since the estimated Reynolds number is #15. This aligns with the classication of "steady laminar inertial ow" for packed beds as outlined by Dixon et al. 55 Diffusion transfer Inter and intra-particle diffusion transfer is a key parameter in chemical systems, as it is inherently connected to the reaction regime. Diffusion transfer within a catalytic bed follows three distinct mechanisms, as described by Krishna: 56 (a) bulk diffusion, (b) Knudsen diffusion, and (c) surface diffusion. Bulk diffusion refers to the transport of species within the empty space of the porous network. It is governed by the relative motion of individual species of a mixture, i.e. molecular collisions, based on their concentration gradients and diffusion coefficients. 56,57 Knudsen diffusion refers to diffusion transfer of species in the proximity of the catalyst surface or in intra-particle space, and is governed by molecule-catalytic surface collisions. 56 Finally, surface diffusion refers to the adsorption and diffusion/migration of molecular species on the catalyst active sites. 56 As all three diffusion mechanisms impact the ow prole and the reaction mechanisms, they should be implemented in the CFD solution.
The porous medium model of Fluent does not consider the existence of solid catalytic particles. 58 Each computational element is considered to be an entirely homogeneous "uid" zone with a momentum sink applied, as a result of the existence of pores, rather than having dedicated uid and solid subzones. Unfortunately, this restriction does not allow the denition of intra-particle space.
Catalytic particles smaller than 0.35 mm, and especially of the size considered here, 0.2 mm, do not suffer from intra-particle diffusion limitations. 13,44 As such, the reaction is governed by the kinetics, where external mass transport is key, 59 so a volumetric diffusion transfer approach, such as the one offered by Fluent, is not unreasonable. For an accurate and realistic ow prole, however, both bulk and Knudsen diffusion transfer must still be considered; it has been observed that Knudsen diffusivity predominates at lower pressures (10 bar), while bulk diffusivity predominates at higher pressures (100 bar). 44 In the considered experimental runs of Park et al., 46 the operating pressure is 50 bar, so both diffusion mechanisms will contribute. Regarding the third diffusion mechanism, surface diffusion, the studied KMs were dened by considering the adsorption of species on active sites. Thus, even though migration of adsorbed species on neighbouring sites is not individually accounted for, surface diffusion mechanisms are already implemented in the reaction kinetics.
While Fluent's code considers bulk diffusivity, 60 its user-interface (UI) does not directly allow the combination of bulk and Knudsen diffusivities. Dixon et al. [61][62][63] carried out a discrete element CFD model which combined the bulk and Knudsen diffusivities in the form of the dusty-gas diffusion correlation, initially introduced by Hite and Jackson. 64 This was possible through the use of Fluent's scalar equation theory. A similar methodology to that of Dixon et al. [61][62][63] is considered here as well. The methodology to implement it into the CFD model is presented in detail in the ESI. †

Chemical reactions
Chemical reactions are applied volumetrically through User-Dened Functions (UDFs), i.e. external user codes, as a source term in the conservation equation of the scalar quantities (see ESI †). This source term is dependent on the reaction rate dened in the considered KMs of Park and Nestler, and follows the general equation form for chemical reactions 60 in Fluent: where r i takes the form of eqn (4)-(6) in the case of Park's KM or of eqn (7) and (8) (9), is applied in units of kg m À3 s À1 , so the reaction rate constants calculated from the respective KMs should be modi-ed accordingly.
The reaction constants follow the form of Arrhenius equations, and their values are presented in the ESI † for both KMs. The equilibrium constants for the reaction rates, K eq,1 , K eq,2 , and K eq,3 , are taken from Graaf and Winkelman. 65 They are also presented in the ESI. † In Fluent's code, it is common practice to dene the carrier gas, argon in this case, as the bulk species. 58 The bulk species mass fraction is always such that the mass balance of species is preserved, meaning that the sum of all species mass fractions doesn't go above unity. For the scalar quantities, however, the conservation law is not ensured, and it must be enforced, as both diffusion and chemical reaction sources are applied directly on the scalar equations. Two mass balance conservation conditions are considered at the end of each iteration, one which forces the sum of all 6 scalars (4 CO 2 , 4 CO , 4 H 2 , 4 H 2 O , 4 MeOH , 4 Ar ) to be equal to unity and a secondary one which forces the reactant and the product scalars to be equal to the initial mass fraction of the feedstock (4 CO 2 /4 CO /4 H 2 ). Aer the conservation law is applied, the species mass fractions are replaced with the values of the scalar quantities, i.e. Y i ¼ 4 i , to complete the iteration.
The experimental runs considered from Park et al. 46 are presented in Table 1. The stoichiometric number, SN, and the carbon oxide ratio, COR, are dened as follows: 53 COR ¼ A selection of cases taken from Park et al. 46 were chosen to focus on. The selected feeds cover a range of SN and COR numbers. Feeds with 0 mol% CO 2 were not considered, as the model of Nestler was derived by excluding these feeds.

Model validation
Park et al. 46 presented their results as CO 2 and CO conversions at the outlet of the reactor. In addition, they identied whether equilibrium was reached for the feed gas compositions or not. This information can be used to validate the two kinetic models. The comparative results, presented as the experimental and predicted CO 2 and CO conversions at the outlet, are shown in Fig. 1. In Tables A3 and A4 of the ESI, † the expected outlet mass fractions of CO 2 and CO, respectively, based on the experimental conversions, are presented, alongside the respective predictions and relative errors of the two KMs. For feed 13, calculating the CO conversion is not applicable, as the initial CO content is 0; the presented conversions are calculated based on the difference between the nal and the initial mass fractions of the species.
Interestingly, both KMs present a very good agreement with the experimental results, for a range of different feed compositions. Overall, Nestler's KM has a smaller error in its CO 2 prediction and a larger error in its CO prediction, compared to Park's KM (Tables A3 and A4, ESI †). Specically, Nestler's KM tends to overestimate the CO conversion. Larger overall errors for both models appear only in feed 25, where the SN increases. It appears that while both models have a similar accuracy when the CO 2 /CO ratio changes, this is not the case for when the SN changes. Still, however, the error for both models is well below 10% in all cases. Nestler's KM was dened for an SN ratio of up to 3.5, 53 so feed 25, with an SN ratio of 5.6, falls outside its denition range. Therefore, the result that Nestler's KM has a smaller error compared to Park's KM in feed 25 further supports the authors' claim that this model can be accurately applied to a wide range of industrially relevant operating conditions.
According to the experimental results of Park et al., 46 feeds 13 and 15 do not reach equilibrium, while the remaining three feeds do. To identify if the models can predict equilibrium as well, the mass fractions of the species along the axial direction of the reactor are investigated, as predicted from the two KMs. These are all presented in the ESI, † in the dedicated results section for each feed. The equilibrium depth varies along the bed as a result of the feed composition and of the considered KM. For feed 13 (Fig. A2 in the ESI †), Park's KM predicts that equilibrium is reached early in the bed, while in Nestler's KM, the species have almost reached equilibrium by the exit of the reactor, as indicated by the upward trends of MeOH and of H 2 O. On the other hand, for feed 15 (Fig. A13, ESI †), the opposite is true; Nestler's KM predicts an early equilibrium depth, while equilibrium is not reached in Park's KM. In most other cases, equilibrium is reached for both models, but Nestler's KM predicts equilibrium slightly earlier in the catalytic bed compared to Park's KM.
It is noteworthy that while the exit compositions predicted from the two KMs are almost identical, the predicted species proles within the reactor differ signicantly. As such, a simple prediction of the exit CO 2 /CO conversions and whether or not equilibrium is reached leads to a valuable comparison. Yet, without noninvasive in situ composition tracking of the experimental setup, it is hard to validate the two proles. The two models reach comparable outlet results by following discrete reaction mechanisms, and the disparities in the CO prediction are not sufficient to prefer one model over the other. For this purpose, the ow and reaction proles within the catalytic bed should be examined and compared for the two KMs, which is an additional benet offered by CFD models. Through this investigation, the predicted favourable mechanistic pathways can be evaluated according to kinetic studies available in the literature. This can evaluate the validity  of the considered reaction mechanisms and identify potential weaknesses in the available models. For this investigation, two different feed compositions were selected, feed 23 and feed 13. The latter is important to consider as well, as pure CO 2 feeds with no initial CO content are important for CCU studies.
Flow prole within the catalytic bed: feed 23 The species mass fractions along a 2D plane surface in the centre of the reactor are presented in Fig. 2 for feed 23. In addition, the rates of reactions R 1 -R 3 , in terms of CO 2 (R 1 and R 2 ) and CO (R 3 ) rate of change, are presented in Fig. 3. Here a positive rate of change is associated with the production of that species, whereas a negative rate of change refers to its consumption. Specically for R 2 , a negative rate is associated with the RWGS pathway (CO 2 consumption). Nestler's KM predicts that the rate of CO 2 hydrogenation and the rate of FWGS reaction are of similar magnitude. The species reach equilibrium early in the bed, with CO almost instantly compensating the amount of CO 2 consumed for MeOH synthesis.
In Park's KM, for CO, its conversion to CO 2 is a more favourable pathway than its direct hydrogenation to MeOH. Surprisingly, the mass fraction of H 2 O follows an inverse relationship to that of CO 2 ; H 2 O reaches a peak value early in the bed, followed by its consumption, along with CO, in the FWGS reaction to supply enough CO 2 for MeOH synthesis. The rate of CO 2 hydrogenation is quite high near the entrance, thus producing MeOH and H 2 O. This leads to a big spike in the maximum H 2 O mass fraction, which is almost triple the respective value predicted from Nestler's KM (Fig. 2). This sudden peak in H 2 O inhibits MeOH synthesis from CO 2 and drives the FWGS reaction, as seen from the relative magnitudes of the rates of change of reactions R 1 and R 2 . The inhibitory effect of H 2 O in the MeOH synthesis reaction is well known, as it has a detrimental effect on the equilibrium. 4,7,23,42 In this case, the outlet mass fractions of the different species predicted from the two models have a much larger variation compared to those of feed 13. Nestler's KM predicts 5% less CO and H 2 O and 3.1% more MeOH at the exit of the reactor, compared to Park's KM, so the inhibitory effect of H 2 O is not as pronounced in Nestler's case. The local temperature variations in the centreline of the reactor, for feed 23, are presented in the ESI (Fig. A19 †). The temperature values correspond to the local variation from the 523 K operating temperature. In addition, the energy released specically from chemical reactions is also shown ( Fig. A20 and A21 †), with positive energy release associated with an energy source or exothermicity, and negative energy release associated with an energy sink or endothermicity. Both KMs predict a similar temperature prole, with a local temperature increase early in the bed, where the magnitudes for the exothermic reactions, MeOH synthesis-R 1 and FWGS-R 2 , are the highest. Interestingly, Nestler's KM predicts a higher local temperature variation prole, despite omitting CO hydrogenation, which is the most exothermic reaction. Park's KM predicts that the magnitude of the CO hydrogenation rate is much smaller compared to the other two reactions early in the bed (Fig. 3), however, due to its signicant exothermicity (À90 kJ mol À1 ), its contribution as an energy source is even larger than that of the FWGS reaction (Fig. A20 in the ESI †). As with all feeds, with the exception of the temperature variation predicted from Nestler's KM for feed 13, the hotspot location is early in the catalytic bed.

Flow proles within the catalytic bed: feed 13
The species mass fractions along a 2D plane surface in the centre of the reactor are presented in Fig. 4 for feed 13, which has no CO in the initial feed. The rates of  reactions R 1 -R 3 , in terms of CO 2 (R 1 and R 2 ) and CO (R 3 ) rate of change, are presented in Fig. 5. Park's KM predicts that CO 2 hydrogenation is more important early in the bed, causing MeOH to reach a peak value before being consumed to produce CO. Specically, CO is produced evenly from RWGS (52.87%) and MeOH decomposition (47.13%), while no CO is consumed in this case to produce MeOH. However, Nestler's KM predicts that initially the RWGS reaction dominates with a higher rate than CO 2 hydrogenation. In this KM, the peak CO concentration is predicted to be early in the bed, followed by its consumption in the FWGS reaction. Unsurprisingly, in both models, at 0 mol% initial CO concentration, production of CO takes priority, as opposed to MeOH production. This conrms the observations in the literature that low CO mixtures are not as effective for MeOH synthesis. 13,39 Kunkes et al. 66 showed, through H/D isotope substitution experiments, that CO is not produced from MeOH decomposition, but instead only through the RWGS reaction. Though in their work, the feed gas composition included both CO and CO 2 with concentrations of 6% and 8%, respectively. In all other feeds where CO is initially present in the feed gas, CO is not produced from MeOH decomposition but only through the RWGS reaction. As such, the possibility that MeOH decomposes to form CO, in the absence of CO, cannot be dismissed. A parametric study to investigate how the initial mass fraction of CO in the feed affects the rate of MeOH synthesis/decomposition for R 3 will be presented in a following section.
The local temperature variation in the centreline of the reactor, for both KMs, is presented in the ESI (Fig. A3 †). The energy released by the chemical reactions is also shown in Fig. A4 and A5 of the ESI. † As predicted from Park's KM, a large amount of energy is released early in the bed, due to the exothermic nature and signicant magnitude of R 1 . As the rates of the RWGS and of MeOH decomposition for CO production, i.e. reverse R 3 , are much smaller compared to that of R 1 , their endothermicity cannot balance the exothermicity of R 1 , causing a temperature increase early in the bed. Deeper into the bed, the reaction rate magnitudes of all three reactions drop, along with the local temperature variation. Nestler's KM predicts a signicantly different local temperature variation prole; since the RWGS reaction, R 2 , is the only reaction which produces CO in Nestler's KM, its magnitude is larger early in the bed compared to that of MeOH synthesis, R 1 . This results in an overall energy sink and a local temperature decrease. When examining the 2D plane surface of the temperature variation, Fig. A3, † uneven local variations were observed in Nestler's case; specically, variations where both the radial and the axial temperature present random uctuations that don't match the species ow prole (Fig. 4) and the reaction rate prole (Fig. 5). Such proles are usually associated with numerical errors, mainly caused by mesh quality. Yet, all other feed compositions for both KMs, including feed 13 with Park's KM, present a uniform temperature variation prole with a smooth transition, which follows their respective species ow proles. To investigate further the reason for this numerical error, and to eliminate the possibility that mesh quality was the reason for the inaccuracies, a mesh independency study was performed for feed 13 as well. Following the same procedure as the mesh independency study presented in the ESI, † a total of 5 different meshes are  investigated, with total mesh element counts of roughly 100k, 300k, 600k (baseline), 900k and 1200k. The temperature variation proles along the 5 different meshes are presented in Fig. 6. In all ve cases, Park's KM presents a smooth temperature transition along the bed. In Nestler's KM, however, the uctuating temperature prole not only persists but also changes across the different mesh cases. It is important to identify the source of the temperature prole difference between the two KMs, as it can be caused by a variety of reasons. One reason could be the existence of reverse ow at the outlet of the reactor, observed in Nestler's KM (Fig. 6), which affects the local temperature variation. A description of what the reverse ow is and how it affects the ow is presented in the mesh independency study in the ESI. † In sum, it is limited only in the nearoutlet region, above the 0.9 axial length of the reactor, and its impact is negligible in the ow prole of the main catalytic body. Thus, it cannot be the reason for the inconsistent temperature prole observed along the catalytic bed. In addition, poor mesh quality can also be excluded as a potential cause, as the mesh quality of the 600k mesh is sufficient to minimize numerical errors (see ESI †). These conclusions allow us to eliminate the possibility that the temperature uctuations are caused by bad mesh quality or incorrect model setup. Their potential cause can thus be limited to the KM itself.
The mass fractions of all species and in all feeds in the outlet are almost identical to their respective experimental values, signifying that the kinetic constants are accurately dened for both KMs. Yet the ow proles of the species, especially of CO and MeOH, differ signicantly between the two KMs (Fig. A31, ESI †). This difference is key, as it originates from the considered reaction mechanisms and causes two distinct reaction equilibrium, as well as local temperature, proles. The only main mechanistic difference between the two KMs, then, is the existence of R 3 . The local temperature variation is dependent on the local balance between endothermic and exothermic reactions; as feed 13 has no initial CO content, the equilibrium of the involved reactions is forced towards its production. In Park's KM, two different reactions are responsible for CO production, RWGS and MeOH decomposition, which both take place simultaneously. In Nestler's KM, however, the exclusion of a direct pathway between CO and MeOH, i.e. R 3 , means that the local species equilibrium (CO 2 , CO and MeOH) is dependent on a single intermediate species, CO 2 . Undoubtedly, the accurate prediction of the correct temperature prole is crucial for large-scale reactors to avoid catalyst sintering, so identifying the error is key for future modeling investigations. This is another instance where CFD can act as a guide for future experimental investigations. One of these two temperature proles is unrealistic, and that can be experimentally conrmed. Following the experimental investigation, a more educated decision can be made on which of the two KMs results in a more realistic species and temperature prole and which mechanisms are key for feeds with no initial CO content.
MeOH production from CO 2 or CO?
Including an expression for CO hydrogenation in the KM can be benecial for investigating questions that have been set in the literature regarding its overall role in the MeOH synthesis process. Specically, it quanties the amount of MeOH synthesised directly from CO. Furthermore, preferable pathways for the conversion of CO into MeOH, i.e. direct hydrogenation or a combination of the FWGS reaction and CO 2 hydrogenation, can be subsequently considered.
In Table 2, a comparison of MeOH synthesis from CO 2 and CO is presented, as predicted from Park's KM. The presented percentage values are estimated volumetrically in the entire reactor. It is important to note that in feed 13, MeOH is volumetrically produced from CO 2 and then is consumed to produce CO. So the negative percentage value for CO and the >100% value for CO 2 are associated with the mechanistic pathways of MeOH.
As can be seen, MeOH is produced predominantly from CO 2 , with a percentage of 62.6%, as opposed to 37.4% of it being synthesised from CO. This agrees with the literature studies that identify CO 2 as the main source of MeOH. 23,34,[40][41][42] In the kinetic study of Grabow and Mavrikakis, 34 for the experimental conditions of T: 499.3 K, P: 29.9 bar, SN: 8.53, and COR: 0.47, they predicted 2/3 of MeOH to be produced from CO 2 and the remaining MeOH to be produced from CO. Their prediction is in excellent agreement with this kinetic model, though their reaction conditions differ from those considered here. As can be seen in Table 2, composition variations (with the exception of feed 13) don't signicantly affect the amount of MeOH produced from CO 2 and from CO. Changing the reaction temperature and pressure, however, will denitely affect the relative balance of the three reaction pathways and possibly alter the amount of MeOH produced from CO. More investigation is required with the produced CFD models to see how MeOH production from CO 2 /CO changes with variations in the operating conditions.
Preferable pathway for CO: RWGS or direct CO hydrogenation?
As already mentioned, the preferable pathway for CO can be investigated with the kinetic model of Park, since CO hydrogenation is included in the rate equation.
To evaluate the preferable pathway for CO, the volumetric CO consumption in the FWGS reaction (R 2 ) and in CO hydrogenation (R 3 ) is presented in Table 3. In addition, the CO consumption in R 2 and R 3 along the centreline of the catalytic bed for feed 23 is presented in Fig. 7. The preferable pathway for CO to MeOH is through conversion to CO 2 , through the FWGS reaction, rather than through direct hydrogenation. The main difference between the two reaction pathways is the large magnitude of the FWGS reaction early in the bed, which, as seen earlier, is mainly driven by the fast CO 2 hydrogenation rate. Later in the bed, the two reaction pathways have a similar contribution to CO consumption.
It has been demonstrated in the literature that CO 2 hydrogenation, in both experiments and in models, is undoubtedly more important in MeOH  production. 40,41,52 The KM of Park agrees with this result, yet it predicts signicant CO consumption for MeOH synthesis as well. In several kinetic models, however, such as those of Bussche and Froment and of Nestler, it is assumed that CO hydrogenation is negligible in the overall mechanisms of MeOH synthesis. The kinetic model of Park predicts otherwise, as CO is consumed in almost equal amounts in both the FWGS reaction and CO hydrogenation. In the kinetic study of Grabow and Mavrikakis 34 it was predicted that only 28% of total CO is consumed by the FWGS reaction, while the rest was consumed by CO hydrogenation. Again, it is difficult to extrapolate these conclusions to the results observed here, due to the vastly different experimental conditions considered.
Other studies predict that the kinetics of the WGS reaction are much faster than those of the hydrogenation reaction 13 but that is not necessarily indicative of the preferable CO consumption pathway, or their relative difference.

Connection between initial CO mass fraction and MeOH synthesis
Specically for feed 13, where there is no CO in the initial feed, it was observed that a signicant amount of the produced MeOH (27.5% in Table 2) is consumed to produce CO. Obviously, this is detrimental to the efficiency of a MeOH synthesis system. The produced CFD model can be used as a parametric analysis tool to investigate the minimum required initial CO mass fraction to negate MeOH decomposition and only drive MeOH synthesis. The KM of Park was used for this parametric study. The volumetric production of MeOH from CO 2 and CO, according to the initial mass fraction of CO in the feed gas, is presented in Table 4. In addition, the CO rate of change for reaction R 3 , presented as the consumption of CO, is shown in Fig. 8. MeOH decomposition to CO occurs up until an approximate initial CO mass fraction of around 0.05, or a 2.5% initial mol fraction. Even at a mass fraction of 0.05, early in the bed, and aer the initial spike in MeOH production in the entrance of the reactor, some MeOH is decomposed to produce CO. When the initial mass fraction of CO is 0.06, or 3.0 mol%, the rate of R 3 primarily drives towards MeOH synthesis rather than towards MeOH decomposition. This demonstrates the benets of cofeeding CO and CO 2 , which have also been observed in the literature for both low CO 13,39 and low CO 2 (ref. 22, 34 and 67) initial content. From the results of the different feeds studied here, it seems that MeOH decomposes to CO only at very low initial CO mass fractions in the feed gas. In all other cases, MeOH synthesis takes place exclusively. This agrees with the results of Kunkes et al. 66 where CO production takes place as a result of the RWGS reaction rather than through MeOH decomposition.

Discussion
One of the most important conclusions that can be reached from this study is the mechanisms for MeOH synthesis from CO 2 . In the paper of Wilkinson et al., 43 it was experimentally observed that early in the bed, CO 2 is consumed to produce MeOH and H 2 O through direct hydrogenation. Deeper in the bed, CO 2 and H 2 O reach a steady level and CO consumption replaces CO 2 consumption. In addition, the production rate of MeOH is faster earlier in the bed than later. 43 All these conclusions agree with the mechanisms predicted from both models, as seen from the species mass fractions along the bed; this indicates that both models provide reasonable results and the nal compositions predicted are close both to each other and to the experimental results. This, however, makes validation of the two models much harder, as this validation requires either a dedicated experimental run or general observations in the literature.
For the latter, these observations are hard to use as validation for the mechanisms, since experimental conditions and setups vary. Grabow and Mavrikakis 34 stated that conclusions regarding the main carbon source for MeOH can only be made for specic conditions and cannot be generalised. The same applies to reaction mechanisms, especially since the errors of the various kinetic models increase outside the operating range for which they were dened. 53 In the current model comparison, some discrepancies could be seen in feed 23, where Nestler's model predicted 3% more MeOH and 5% more CO and H 2 O compared to Park's KM. This result could not be validated with the experimental results of Park et al. 46 Actually, this result was fundamentally related to the preferable reaction mechanisms and whether H 2 O inhibition takes place or not. Conclusively, being able to experimentally identify the physical trends of the different involved species is as important as the actual species values under distinct experimental conditions. The experimental procedure followed by Wilkinson et al. 43 is a great example of the former.
Dedicated experimental runs to dene the reaction mechanisms, given in a form that can be usable by various kinetic investigation methodologies, i.e. DFT, kinetic Monte Carlo and CFD, will be crucial for reaching a common understanding. Especially for papers where experimental results are used to establish new kinetic models, thorough investigations and results are required. In addition, combining our available investigative methods can lead to invaluable information regarding the MeOH synthesis process. As an example, Higham et al. 37 predicted that CO 2 could bond dissociatively on Cu(110) and Cu(100) facets, thus initialising CO hydrogenation on these sites. Yet the adsorption of CO feedstock on these facets was not considered. Park's KM considers that CO 2 and CO adsorb non-competitively on Cu sites while Nestler's KM considers competitive adsorption of CO 2 and CO. The role of these adsorption sites cannot be identied with the CFD model but DFT models are better suited for such investigations. Following the identication of the accurate micro-scale mechanisms, more accurate KMs can then be experimentally established, and their validation on both the local and the reactor scale can be veried with CFD tools. Here, the role of CFD was to investigate the kinetic mechanisms of MeOH synthesis from CO 2 / CO/H 2 throughout the catalytic bed. The produced results can now be used to guide the experimental setup, which in turn can investigate the gas compositions at different reactor lengths. This can help validate the two KMs better and increase our understanding of the overall kinetics. Aer further validation, the CFD model can also be used as a predictive tool for the identication of temperature hotspots or for process optimisation, e.g. if the equilibrium depth can be predicted beforehand, this can help minimise the total catalyst load required. In general, collective and combined research with the available investigative methods has the potential to accelerate scientic understanding.
This leaves the question of which of the two KMs is more accurate and produces more reasonable results. Unfortunately, the answer is not straightforward. In their paper, Slotboom et al. 68 compared several kinetic models, including those of Graaf and of Bussche and Froment, to determine their ability to predict experimental results. They separated the kinetic models into four groups, according to the mechanisms the KMs used to predict MeOH synthesis. They concluded that it is not the inclusion of CO hydrogenation that is key to predicting MeOH synthesis, but rather considering at least two distinct active sites, i.e. adsorption of CO on catalytic sites should also be considered. The considered KMs fundamentally consider three (Park's KM) and two (Nestler's KM) active sites, which allows both models to give accurate predictions. One of the criteria set by Slotboom et al. 68 for choosing between models is the number of parameters necessary to dene that model. The model of Park requires 10 parameters to be experimentally dened, while that of Nestler requires only 7. This makes the model of Nestler easier to use and implement.
In reality, though, choosing between these two models depends on the operating conditions and the expected investigation outcomes. The errors of Nestler's KM decrease as the COR ratio increases, so for feed compositions with high CO 2 content, Nestler's KM is expected to give a more accurate prediction, especially for the CO 2 species. In addition, Nestler's KM presented a smaller error in its prediction for the high SN ratio case, compared to Park's KM. As this model was dened to be applicable to a wide industrially relevant operating range, this benet can be a great asset when deciding the most suitable model. Finally, it is easier to adapt to different experimental setups, as it only requires 7 parameters to be dened. However, it is important to consider that for feeds with no initial content of CO, inaccuracies in the predicted temperature prole can occur. The model of Park can more accurately predict feeds with low COR number. In addition, the inclusion of CO hydrogenation in the KM provides a more comprehensive mechanistic understanding of the effects, the interactions, the inhibitions, and the overall equilibrium of the involved reactions and species. In addition, as this model considers three active sites, from a catalyst design/optimization point of view, it can better guide experimental efforts to optimize the catalyst with the goal of promoting a different mechanism, e.g. to enhance the rate of CO 2 or CO hydrogenation. However, its requirement for 10 parameters makes this model's adaptation to different experimental conditions more challenging.

Conclusions
Methanol synthesis from a mixture of CO 2 /CO/H 2 has been widely studied for its kinetics and reaction equilibria. In this paper, two of the available kinetic models in the literature, those of Park et al. 46 and Nestler et al., 53 are compared using a CFD simulation. The kinetic models are validated, based on previously published experimental results, in terms of their mechanistic predictions within the catalytic bed. As has also been observed in the literature, CO 2 is converted early in the bed to MeOH but reaches a steady state soon aer. CO consumption then replaces CO 2 consumption, as it is converted to MeOH either through the FWGS reaction and subsequent (CO 2 ) hydrogenation or through direct hydrogenation. The inclusion of CO hydrogenation in the kinetic model of Park offers more detailed mechanistic insights in terms of the interactions, the inhibitions, and the equilibrium of the involved reactions. However, Nestler's kinetic model offers a broader applicability, as it was dened for a wide range of industrially relevant operating conditions.
The benets of applying a validated CFD model as an investigation method were demonstrated. The presented CFD model offers a realistic representation of the spatial proles of local species concentrations, temperature variations, and reaction rate magnitudes. Data such as this are hard to acquire experimentally, as non-invasive and in situ experimental setups have limited applicability. More importantly, the CFD model can detect internal inconsistencies in the applied kinetic models or key mechanisms for the MeOH synthesis process at low initial CO content, thus paving the way for further dedicated and educated experimental investigations. In combination with these experimental runs, it can be used to increase our understanding of the methanol synthesis process, and in turn, act as a guide for catalyst design and reactor engineering. Especially for large scale chemical reactors, where knowledge of spatial heat and mass transfer is crucial to avoid runaway conditions, accurate CFD models can play a decisive role.

Conflicts of interest
There are no conicts of interest to declare.