Parametric Analysis and Optimization of Vanillin Hydrodeoxygenation Over a Sulfided Ni-Mo/δ-Al2O3 Catalyst Under Continuous-Flow Conditions

A fundamental understanding of the process parameters affecting the catalytic hydrodeoxygenation (HDO) of bio-oils is of significance for enabling further progression and improvement of industrial biofuel upgrading methods. Herein, a novel demonstration and evaluation of the effect of temperature, pressure, and weight hourly space velocity in the continuous HDO of vanillin to cresol over a Ni-Mo/δ-Al2O3 catalyst are presented. Response surface methodology was used as a statistical experimental design method, and the application of central composite design enabled the generation of a statistically significant simulation model and a true optimization parametric study. The distribution of Ni and Mo on δ-Al2O3 was confirmed using scanning electron microscopy with energy dispersive X-ray analysis (SEM-EDX). No gradients with EDX mapping could be identified, and the elemental analysis showed well-dispersion of the metals. The mesoporous character of the catalyst-support system was unraveled using N2 physisorption. Experiments were conducted within the parametric range of 250–350 °C, 3–9 bar, and 15–35 h−1. Both temperature and pressure were found to have statistically significant linear and quadratic effects on the selectivity for cresol. The parametric interaction of temperature with pressure and space velocity also had a significant effect on the resulting response. The optimal temperature range becomes more critical at lower space velocities. Optimal selectivity for cresol was established at 314 °C, 5 bar, and 35 h−1. The fitting quality of the generated regression model was statistically confirmed and experimentally validated to describe the specified HDO process within the 95% two-sided confidence interval.


Introduction
The demand for renewable and environmentally friendly alternatives to conventional fossil-derived fuels is rapidly increasing. The use of fossil fuels has an established negative environmental impact. Numerous reports confirm and derive, for example, the increase of greenhouse gases in the atmosphere and its associated negative effect on nature and on living life from the usage of fossil fuels. It is stressed that further commitments and measurements in regard to developing alternative renewable fuels remain vital [1][2][3]. Today's commercially available bio-liquid fuels, such as biodiesel and bio-alcohols, are processed from established technologies and are attractive fossil-substituted fuels. However, these are processed from the so-called first-generation feedstocks, i.e., edible commodities (e.g., sugars, starch, and vegetable oils), and contribute to other associated negative environmental impacts [4,5].
A promising second-generation (non-edible) resource for producing biofuels is lignin [6]. Lignin is a complex oxygenated aromatic polymer comprising almost a third of lignocellulosic biomass and is one of the most abundant biopolymers worldwide [6][7][8]. Lignin is produced in large amounts as a by-product in the pulp and paper industry [6,8,9]. Due to its heterogeneity and complex multifunctional structure, its current use is mainly restricted to combustion and on-site energy recovery [8,10]. However, by enabling the utilization of lignin's energy-rich structured monomers, its environmental and economic value would increase and possibly enable its use as an alternative bulk biomass feedstock for bio-based fuels [6,7,9,11].
Isolated biomass feedstocks such as lignin need to be upgraded via appropriate transformation methods prior to being adaptable for subsequential utilization as, e.g., biofuels. Several techniques have been established to convert biomass feedstocks into bio-oils, including fast pyrolysis, gasification, hydrogenolysis, oxidation, and hydrolysis [12]. Fast pyrolysis, which is a thermal degradation process performed under anaerobic conditions, is the most recognized and used technique. It is considered beneficial due to being flexible in terms of feedstocks and for being cost effective [12,13].
However, a general challenging aspect of using biomass as feedstock in comparison to fossil-derived hydrocarbons is its difference in composition. The major compositional difference relates to the extensively higher content of thermodynamic unstable and reactive oxygen-containing molecules in bio-oils derived from biomass compared to conventional oils [1,3,4,14,15]. Thus, bio-oils need to be stabilized before further utilization [1,2,4,6,14,15].
Catalytic hydrodeoxygenation (HDO) is an established hydroprocessing method for bio-oil upgrading and is considered the most promising oxygen removal technology in this respect [1,2,14,16]. HDO is adapted from industrial heteroatom removal hydrotreatment processes associated with fossil-derived fuels, e.g., hydrodesulfurization and hydrodenitrogenation [3,14,15]. However, major challenges are related to the oxygen removal hydrotreatment process [3,15]. The HDO reaction is exothermic, and the extensive amount of highly reactive oxygenated containing molecules in bio-oils results in a significant risk of an accelerated deactivation of the catalyst due to, e.g., associated polymerization and resulting coke formation [4,14,15]. Furthermore, the large size of bio-oil-containing molecules limits the diffusion of the reactants and products in narrow pores and channels, which further complicates bio-oil upgrading [4,15,17,18].
The catalyst-support system used in HDO is decisive for the resulting converted liquid-fuel quality. It is a complex area with several synergistic aspects. The catalytic performance and the reaction pathways are generally dependent on metal dispersion and agglomeration, chemical properties of active sites (e.g., acidity/basicity), and on surface properties (e.g., pore volume and surface area) [5,6]. Different reaction pathways for the HDO of bio-oils are reported in the literature, e.g., direct hydrodeoxygenation, decarbonylation and decarboxylation, and are dependent on the feed and the catalyst [5,19]. Comprehensive research in this regard is reported, and a variety of catalysts have been reported to be active in the HDO process, e.g., transition metals, noble metals, and phosphides [2,3,6,15,20]. The most reported and used catalysts in hydrotreatment processes are sulfided Mo promoted with Ni or Co on γ-Al 2 O 3 support (Ni-Mo/γ-Al 2 O 3 or Co-Mo/γ-Al 2 O 3 ) [3,4].
Efficient and controlled operational conditions in the HDO process are challenging to achieve due to the wide variety of oxygenated multifunctional groups in bio-oils [4,6,15]. The reactivity differs between functional groups, and the resulting suitable temperature and pressure vary between groups and certain compounds [4]. Hence, complete hydrogenation, polymerization, coke formation, and diversity of resulting products are significant problematic aspects to consider in sufficient HDO processes [3,15]. To enable more efficient operational conditions for bio-oil upgrading processes, the hydrotreating is commonly performed by co-processing with petroleum fractions or in a two-step hydrotreatment process. The two-step hydrotreating process includes an initial stabilization step, partial deoxygenation, performed under milder conditions and a subsequential step, under more severe conditions, aiming at complete deoxygenation [2-4, 14, 15, 21]. Improved methods for bio-oil upgrading are required, and a fundamental understanding of reaction pathways and process affecting parameters are of significance [2,15].
The valorization of lignin is challenging, and model compounds are commonly used to gain a deeper understanding of the reaction pathway and parameters affecting the performance [6,8,14]. Several studies on the HDO reactions, substantially performed in batch-reactors, of lignin-model compounds, have been reported elsewhere [6,15]. However, to develop improved methods for bio-oil upgrading HDO processes in the industry, a more fundamental understanding of reaction pathways under continuous-flow conditions is important [15]. General affecting factors such as possible mass-transfer gradients and reaction rate-limiting steps are different in continuous-flow packed-bed reactors compared to batch reactors [22,23]. Studies performed under continuous-flow conditions are also mainly limited to lignin-model compounds with one or two oxygen-containing functional groups, e.g., cresol and guaiacol, both being phenolic compounds, which restricts the understanding of HDO process pathways [1,8,[14][15][16].
In the present study, vanillin was chosen as a representative oxygenated model compound for the continuous catalytic HDO to cresol over a Ni-Mo/δ-Al 2 O 3 catalyst. Vanillin is an aromatic compound containing several oxygen groups (ether, hydroxyl, and aldehyde groups) that are present in bio-oil, which makes it a suitable model compound in that 1 3 respect [1,16]. Central composite design (CCD), as part of response surface methodology (RSM), was used as a statistical design of experiment method to demonstrate the optimal operating conditions in respect of temperature, pressure, and space velocity. The continuous catalytic HDO system is a multivariable system, and by using CCD, the number of required experiments to give a comprehensive understanding of the affecting parameters is reduced.

Catalyst Synthesis
The δ-alumina-supported nickel-molybdenum catalyst (Ni-Mo/δ-Al 2 O 3 ) was prepared by incipient wetness impregnation. δ-alumina extrudates (Sasol Germany GmbH, Germany) were crushed and sieved into a particle size range of 1-2 mm. The metals were impregnated separately in sequential steps. An aqueous solution of molybdenum (ammonium molybdate tetrahydrate, ACS Reagent, Sigma-Aldrich) was impregnated initially, followed by drying at 120 °C for 4 h and calcination at 500 °C for additional 4 h (ramp 120 °C/h in both steps). An acidic aqueous solution of nickel (nickel(II) nitrate hexahydrate, Sigma-Aldrich) was impregnated sequentially, using the same method and the same drying and calcination program. Citric acid monohydrate was used to prepare the acidic solution, using a 0.7 molar ratio with respect to nickel.

Catalyst Characterization
Key characteristics, in terms of metal distribution and surface properties, of the synthesized Ni-Mo/δ-Al 2 O 3 catalysts were examined by using scanning electron microscopy (SEM), and by N 2 physisorption analysis, respectively.

Scanning Electron Microscopy
High-resolution images and elemental analysis of the synthesized catalyst were performed using a field emission scanning electron microscope, JEOL JSM-6700F (JEOL, Akishima, Tokyo, Japan) instrument equipped with an energy dispersive X-ray system (EDX). The associated field emission gun was a X-Max model (Oxford Instruments, Abingdon, UK). An approximately 40 nm thick layer of carbon was deposited on the catalyst prior to analysis to enhance the conductivity of the surface.

HDO Activity Test and Procedure
A parametric investigation of the catalytic HDO of vanillin to cresol over Ni-Mo/δ-Al 2 O 3 was performed in a continuous-flow fixed-bed reactor. The process system design, schematically illustrated in Fig. 1, included the inlet of feed, Fig. 1 Schematic illustration of the continuous-flow reactor setup used in the HDO experiments. MFC is the abbreviation for mass flow controller, whereas T1 is a temperature controller and P1 and P2 are pressure regulators pre-heater, reactor, the condensation step, separation of gas and liquid, and a liquid collecting vessel. The liquid reactants were fed using a HPLC-pump, LC-10ADvp model (Shimadzu, Kyoto, Japan), whereas the gas was fed by using a mass flow controller (MFC) (Bronkhorst® High-Tech BV, Ruurlo, the Netherlands) and mixed in the pre-heater. To control the pressure (gauge) and temperature in the system, a pressure controller and temperature controller were used. The reactor was made of a stainless-steel tube with an inner diameter of 17 mm. To reduce possible mass and temperature gradients in the reactor bed, the catalyst was mixed with the inert material α-alumina (10-20 mesh) with a volumetric ratio of 1:1. An equivalent layer of the same inert material was also added both below and under the catalyst bed to further reduce possible gradients. The amount of catalyst used was 1 g. Technical xylene (VWR Chemicals, ≥ 98%) was used as a solvent, and a liquid feed of 7400 ppm vanillin (TCI Chemicals, > 98%) was used throughout the experiments. The gas, 1% hydrogen sulfide (H 2 S) in hydrogen, was fed with a molar ratio of 10:1 with respect to vanillin. To reduce the risk of product contaminations from previous runs, the liquid collecting vessel was emptied after 90 min, whereafter the produced liquid during the subsequent 30 min was collected for product analysis. Before performance tests, the catalyst was sulfided in the reactor at 400 °C for 4 h (ramp rate to temperature from ambient 120 °C/h) at 8 bar using 1% H 2 S in H 2 (WHSV 400 h −1 ).

Product Analysis
The produced liquid was analyzed using gas chromatography, 456-GC (Scion Instruments, Goes, The Netherlands), coupled with a flammable ionization detector (FID) and a Rtx® DHA-50 column, 50 m × 0.2 mm (Restek, Centre County, PA, U.S). The conversion (Conv.%) of vanillin was calculated according to Eq. (1), whereas the selectivity (Sel.%) to cresol was calculated according to Eq. (2). The instrument was calibrated with standard chemicals purchased from Sigma-Aldrich.

Design of Experiments
The multivariable parametric effect on the performed catalytic HDO process was investigated by using the response surface methodology (RSM) and the CCD. RSM includes statistically designed experiments, mathematical modeling, resulting response prediction models, and validation for the modeling accuracy of the system [26][27][28][29]. The CCD allows for an accurate evaluation of parametric effect and enables a true optimization analysis of the system [29,30]. The regression modeling and analytical evaluation are based on the coding of the variables. By introducing coded dimensionless variables, + 1 to − 1, the resulting importance, i.e., the effect, that each factor has on the response, can be directly relatively compared and determined [26,28,29]. The high level of each specific variable is donated + 1, and in contrast, the low level -1. The addition of a central point (coded as 0) and surface axial points (denoted ± α) in the design of experiments, as done in CCD, enables the detection of possible curvature in responses [26,29,31]. Using axial points outside of the surface level enables the detection of possible optimal conditions outside of the set parametric range. The conventional factorial experiment method, one-factor-at-atime, limits the process evaluation to each of the affecting factors and not the combined affecting interactions [31]. An assumed optimization point for the chosen response is restricted to only the actual performed experiments data. By using CCD, simultaneously, parametric change is conducted, and an optimization modeling is maintained over the whole parametric range [26]. By modeling the response, statisticbased prediction regression models can be established and further used for evaluating the outcome responses for different parametric changes. The experimental points in CCD are schematically illustrated in Fig. 2.
The independent process parametric variables investigated in this study were temperature, pressure, and weight hourly space velocity (WHSV). The parametric intervals of  Table 1.

Modeling and Optimization
Analysis of the experimental response data was conducted using the Design-Expert® software, version 13, Stat-Ease, Inc., Minneapolis, MN, USA, and is statistically based on the analysis of variance (ANOVA) methodology. ANOVA is a statistical tool for evaluating regression models and the statistical significance of the model fit (e.g., F-and p-values). The software was used to identify parametric statistical significance, regression modeling and optimization, and to graphically plot the resulting data. The optimization analysis was conducted with respect to the responsive selectivity for cresol.
In order to minimize the negative environmental impact of the process, such as the associated increased need for energy at higher temperatures and the amount of feed at increased feeding rates, the predicted regression equation was proposed in respect of minimizing the parametric importance, while still maximizing the response of produced amount of cresol. The experimental statistical significance validation of the resulting regression model for the system and production of cresol, was examined by performing experiments after re-loading of catalyst at proposed optimized conditions three times.

Predictive Modeling and Regression Analysis
The conversion of vanillin is consistently high for all the performed experiments, indicating the high activity of the catalysts during the reaction. The generated experimental CCD test matrix and the resulting product analysis are given in Table 2.
Selectivity data to cresol was used for parametric analysis and optimization. An initial cox-box transformation of the experimental data was suggested to square rot. A cox-box transformation is an accepted tool that allows for normality of the data and reduces risks for model factor bias [30,32]. The resulting generated coded quadratic regression model is given in Eq. (3), where A, B, and C represent the coded variables in terms of the process parameters temperature, pressure, and WHSV, respectively.
The variance inflation factor (VIF) gives an estimate for the linear independent effect on the response of a specific factor in a regression model with multiple factor terms. Orthogonality results in VIF equal to 1, and a linear dependence of the specific factor on the response is confirmed [33,34]. Linear dependence of model factors A, B, C, AB, and AC was confirmed (VIF = 1). VIF values greater than 1,  imply multi-collinearity dependency of other model factors for that specific factor. However, VIF values below 5 are generally tolerable without the need for further correction [33,35]. The indication of multi-collinearity existence for the model factors A 2 and B 2 (VIF = 1.05) was hence considered low and negligible, and the effect of each specific model factor, in sequential evaluation, was considered to have an accurate linear dependency on the response. Furthermore, coding of the variables enables, e.g., dimensionless factor comparison and a direct understanding of the relative importance, i.e., the effect of the factors on the response [26,28]. The resulting factorial effect on the responsive outcome is graphically presented in the Pareto plot in Fig. 3. The relative factorial importance of the response can also be distinguished from the generated model. The factorial effect on the response represents by a two-unit change per factor, whilst the corresponding change for the coefficient is represented by a one-unit change [31].
The Pareto plot, Fig. 3, graphically presents the effect of each factor on the response separately and with descending relatively importance from left to right. The relative percentual effect from each factor affecting the total response is also presented. The quadratic interference of temperature, model term A 2 , has the most influence on a responsive change with respect to the other factors, and the effect corresponds to 30% of the total change in response. A change of model term C (WHSV) does not have a high effect on the response, relative to the other factors. However, the combination of temperature and WHSV (AC) influences the response.
The statistical ANOVA significance evaluation of the fit quality of the generated model (Eq. 3) and model terms are presented in Table 3. Table 3 presents the resulting calculated sum of squares (SS), mean of squares (MS), degrees of freedom (DF), F-values, and p-values for the analyzed generated model. The SS and MS demonstrate if the variance of the source has a measurable effect, whereas the F-value and p-value indicate if the effect of variance is statistically significant. The presented F-value is the observed value, which, compared to tabulated critical reference distribution values for a specific confidence interval, demonstrates the chance of the effect being due to the generated model and  model terms rather than test system noise (source of error). The calculated p-value indicated the probability that the generated model and included model terms have a statistically significant effect within a specific confidence interval [30,31]. 95% confidence interval was used in this evaluation. The high F-value of 84.71 implies the statistical significance of the generated model and a 0.01% chance that the model fit is due to noise. p-values less than 0.05 indicate statistical significance [26][27][28][29]. Hence, the resulting low p-value for the model (< 0.0001) also statistically implies that the model quality is of significance [27,34]. In accordance, the lack of fit of the model is stated as not significant (p-value > > 0.05). Furthermore, significant model terms, in this case, are A, B, AB, AC, A 2 , and B 2 . The resulting accuracy statistics of the generated regression model to the experimental data are presented in Table 4.
The adequacy precision measures the signal-to-noise ratio and indicates the suitability of using the model for process guidance [27,29]. Adequacy precision values greater than 4 are desirable and indicate a suitable predictive model for the real system [27,29,34]. The resulting value of 26.03, in accordance with previous ANOVA model fit quality data, implies that the model can be accurately used for parametric changes and understanding of the experimental process outcome. The coefficient of determination, R 2 , describes the part of total variation that is given by the model and gives an estimate of how accurate the model describes the data [27,29]. R 2 value close to one is desirable [27,30]. The R 2 value for this model was 99%, indicating that the data variation is desirably described by the model. Regression models, perfectly fitted (R 2 = 1), are however not necessarily representative and accurate. The compatibility of the generated model can always be adjusted and improved to fit all data points by adding higher ordered model terms. A perfectly fitted model can be a result of the inclusion of data noise, and the predictive system model can give an inaccurate understanding of the system and possible resulting data. In this respect, the comparison of the R 2 to the given adjusted and predicted R 2 is important. Contrary to the R 2 being improved by the addition of model terms, the adjusted R 2 accounts for the adding value of the additional modeling terms. Hence, the addition of not significant terms in the model results in a decrease in the adjusted R 2 value [30,31]. The adjacent high resulting values of the R 2 and the adjusted R 2 , 99% and 97%, respectively, imply a model including sufficient value-adding terms. Furthermore, the predicted R 2 gives an estimate of the accuracy of the model to observations. In accordance with the adjusted R 2 , the predicted R 2 decreases in relation to the increase of included noise affecting data in the regression model [27]. Hence, the accuracy of the model to the system is additionally confirmed with the predicted R 2 of 92%. In addition, the low, resulting standard deviation of 0.23 also indicates that the model gives a predicted value adjacent to the real responsive value.
The statistical fit quality of the generated regression model was improved by reducing the statistically not significant effecting model product terms (BC, C 2 ). The independent model term of C (WHSV) was, however, not removed. WHSV is an inevitable process parameter and was hence considered to be an important specified model term. Additional higher-order product terms were also considered unnecessary for possible improvement of the model fit.
Furthermore, the resulting data for the experimental run, performed at a pressure above the specified parametric pressure interval ( Table 2, run 5), was considered not reliable due to problematic system limitations. At this condition, the liquid was condensed prior to the actual condensation step, causing questionable accuracy of the remaining product composition. The results from this experiment were, hence, not used for further regression modeling and parametric evaluation. To ensure the desirable orthogonality of the coded matrix, the result from the experiment performed, in contrast, with a pressure below the specified parametric pressure interval ( Table 2, run 6), was not either used for further evaluation. The resulting VIF values of the model factors confirmed the preservation of the test matrix orthogonality. An accurate evaluation of the linear association effect of the specific factor on the response was preserved.
A final observational analysis of the quality of the model is shown in Fig. 4.
The graphically presented data in Fig. 4 further complies with the accuracy of the generated model and its predictive use for the system. In Fig. 4a, the actual experimental data is plotted versus the resulting predicted data when using the model. As can be seen, a strong correlation between the actual and the predicted data was observed with minor deviations from the 45° line. This linear relation indicates that experimental data can be accurately explained by the regression model. Furthermore, in Fig. 4b, the externally studentized residuals are plotted for each experimental data point (run number) and give a graphical understanding of the existence of potential outliers. The externally studentized residuals give an estimate of the relative error, in standard deviation units, between data points [30]. It presents the ratio between the residual and estimated standard deviation

Parametric Analysis and Optimization
Both temperature and pressure have a separately linear and quadratic significant effect on the responsive selectivity for cresol. In addition, the combined parametric interaction effect of temperature and pressure, temperature, and WHSV, also have a statistically significant effect. However, the independent parametric variable of WHSV (C) does not have a significant effect (F-value < 4 and p-value > 0.05) on the selectivity for cresol in this system. The resulting corresponding response analysis, in terms of selectivity for cresol, of using the regression model (Eq. 3) over the process parametric interval is graphically presented in Figs. 5, 6.  Figure 5 presents contour plots of the resulting combined process parametric effects in terms of temperature versus pressure and WHSV, respectively, within the experimentally specified parametric intervals on the selectivity for cresol. Whereas Fig. 6 gives three-dimensional representative response surface graphs. It can, graphically, clearly be distinguished that temperature and pressure have a curvature influence on the resulting outcome and that optimum conditions for the system in respect of selectivity of cresol exist. It can be distinguished that pressure in the lower range, approximately 3-6 bar, is desirable and combined with a temperature in the approximate range of 280-320 °C results in the highest selectivity for cresol. The WHSV does not have a curvature influence on the response. However, at lower WHSV, the temperature range is more critical in respect of the highest possible selectivity for cresol.
Proposed hydrodeoxygenation reaction pathway of vanillin to cresol presented in literature suggest initial hydrogenation of the aldehyde group followed by hydrogenolysis (forming creosol via vanillyl acid) [1,36], and subsequential demethylation with simultaneous hydrogenation and hydrogenolysis (forming p-cresol) [37,38]. The selectivity towards cresol is a synergistic matter, being dependent on, e.g., the catalyst and the process conditions. The catalyst needs to possess suitable properties for the mechanism, e.g., possess acid sites [5]. A certain process temperature is also required for the catalyst to be active. Acid sites on the catalysts and/ or increased temperatures promote cleavage of the C-C bonds [5], and the initial formation of guaiacol instead of creosol is possible. From guaiacol it has been reported that o-cresol can be formed, in addition to p-cresol, by hydrogenation, demethylation and methylation [19,37]. These reactions have been observed to be prone on alumina support compared to demethoxylation. Methylation in the m-position is reported not to be favored [37]. Furthermore, increased pressures of hydrogen promote more efficient hydrogenation and hydrogenolysis [39]. However, increased hydrogen pressure promotes saturation of carbons. High residence times promotes deoxygenation and lower flow rates are recognized as favorable [5].
The presented demonstration of the influence of the process parameters and their combined effect for this system is novel. The results reflect the synergistic affecting factors with respect to temperature and pressure, both having quadratic dependency. The linear dependency of the flow rate in this system could be due the experiments being performed at a high flow rate range in which the flow change was too narrow to make a significant difference in mass transfer.
The statistically validated regression model (see Sect. 3.1 and Eq. 3) was used to find optimum conditions for the highest responsive selectivity for cresol in the continuous catalytic HDO process system used in this study. The final optimum conditions were established after consideration of the parametric environmental impact and still maximizing the responsive outcome. From an environmentally sustainable perspective, less severe process parameters are generally less energy-consuming and hence argued as desirable. The predicted optimal selectivity for cresol (33.4%) in the specified system was established to be at 314 °C, 5 bar, and 35 h −1 .

Experimental Verification of the Regression Model
The resulting optimal conditions for maximum selectivity for cresol (314 °C, 5 bar, and 35 h −1 ) were used for experimental verification of the generated regression model (Eq. 3).

Fig. 6
Three-dimensional representative response surface graphs of the resulting combined process parametric effects on the selectivity for cresol, using the regression model, in terms of a temperature versus pressure b temperature versus WHSV

3
The reactor was re-loaded with a new catalyst bed prior to the experiments, according to the procedure described in Sect. 2.3. The experiments were repeated three times. The resulting selectivity for cresol for the experiments was 30.5%, 31.4%, and 32.7%. All experiments resulted in a conversion of vanillin > 99%. The resulting standard deviation of the selectivity was estimated to be 2.65 within a 95% confidence interval, implying experimental verification of the generated regression model. It also indicates equivalent temperature and mass transportation gradients in the system. The experimental verification results confirm the reliability of the specific process system design used in this study. The regression model can be accurately used as guidance for process parametric changes with respect to temperature, pressure and WHSV in the specified parametric interval.

Catalyst Morphology and Texture
The surface morphology and metal distribution of the synthesized Ni-Mo/δ-Al 2 O 3 catalyst used in the HDO experiments of vanillin, were examined by SEM-EDX and N 2 physisorption. The resulting images of the SEM-EDX analysis of the metal distribution are presented in Fig. 7.
Dispersion of Ni and Mo on the alumina-supported synthesized catalyst was confirmed with SEM-EDX. No metal distribution gradients could be identified with elemental mapping on a surface cut of a particle, confirming a good dispersity of the metals (Fig. 7). Four different positions on the outer surface and the surface cut of a particle were used for quantitative elemental analysis, shown in Fig. 7. The average atomic percentage on the outer surface and the surface cut was estimated from the four investigated positions to be 21 at% Ni and 79 at% Mo, and 38 at% Ni and 62 at% Mo respectively, giving average metal ratios of Ni:Mo to be 0.3 (SD 0.050) and 0.6 (SD 0.026). The differences between the estimated average atomic percentages at the outer surface and the surface cut can be related to either higher Mo impregnation on the surface and/or higher impregnation in the sequential step of Ni in the pores, both giving resulting higher ratio of Ni:Mo in the pores compared to the surface. The low standard deviation values and the resulting EDX mapping do, however, propose uniform impregnation of respectively surface, indicating a well-dispersion of Ni and Mo on the particle.
The specific surface area of the synthesized catalyst was calculated from the resulting N 2 adsorption-desorption isotherms using the BET method and was estimated to be 101.1 m 2 /g catalyst. Using the BJH method, the pore volume and the average pore size were estimated to be 0.6 cm 3 /g and 22.8 nm, respectively, suggesting a mesoporous character of the catalyst support, consistent with the equivalent catalystsupport system reported in the literature [40].

Conclusion
A novel optimization study of the catalytic HDO of vanillin to cresol under continuous-flow conditions has been described. By using the CCD methodology for the experimental design combined with ANOVA for statistical regression, the multivariable parametric effects of temperature, pressure, and WHSV of the process could be accurately determined. Temperature and pressure had statistically significant linear and quadratic effects on the responsive Fig. 7 SEM images, (a, b), and EDX elemental distribution, c-f of the Ni-Mo/δ-Al 2 O 3 catalyst. SEM images show positions for quantitative elemental spectrum analysis. a The outer surface of a particle. b Cross section surface of a particle. c Aluminum. d Oxygen. e Molybdenum. f Nickel selectivity for cresol. The effects of parametric interactions of temperature and pressure and temperature and WHSV on the response have also been concluded to be statistically significant. The process temperature has been shown to be more critical at lower WHSV. Optimal operating conditions for the HDO of vanillin to cresol exist in the studied parametric range, and optimal selectivity has been established at 314 °C, 5 bar, and 35 h −1 . The quality of the generated regression model is statistically significant and can be accurately used in the determination of possible responsive outcomes for different parametric changes, including temperature, pressure, and WHSV. Good dispersion of Ni and Mo on the surface of the catalyst-support system has been confirmed with SEM-EDX. Evaluation of the resulting adsorption-desorption isotherms by the BET and BJH methods suggested a mesoporous character of the catalyst. The findings of this study can be used to gain a deeper understanding of vital affecting process parameters in the continuous HDO of bio-oil-derived oxygenated compounds. Work is ongoing on a further fundamental understanding of the catalyst system to develop improved methods for bio-oil upgrading in industrial contexts.