Optimization Study of Biomass Hydrogenation to Ethylene Glycol Using Response Surface Methodology

: Statistical-based study using response surface methodology (RSM) was conducted to study the e ﬀ ects of process parameters towards biomass hydrogenation. Using Malaysian oil palm empty fruit bunches (EFB) ﬁbres as feedstock, the central composite design (CCD) technique was employed and 18 runs were generated by CCD when four parameters (mass ratio of binary catalyst, hydrogen pressure, temperature and mass ratio of catalyst to feedstock) were varied with two center points to determine the e ﬀ ects of process parameters and eventually to get optimum ethylene glycol (EG) yield. RSM with quadratic function was generated for biomass hydrogenation, indicating all factors except temperature, were important in determining EG yield. Analysis of variance (ANOVA) showed a high coe ﬃ cient of determination (R 2 ) value of > 0.98, ensuring a satisfactory prediction of the quadratic model with experimental data. The quadratic model suggested the optimum EG yield should be > 25 wt.% and the EG yield results were successfully reproduced in the laboratory. and 180 g of 0.438 mol L − 1 of NaOH solution at at for at Upon


Introduction
Using catalysts and agricultural lignocellulosic waste as feedstock, bio-based ethylene glycol (EG) could be produced. A research team in China filed a patent in 2008 indicating ethylene glycol could be produced from cellulose under a catalytic hydrogenation process [1]. Using cellulose, glucose or biomass as feedstock, extensive work has been carried out by researchers around the world who have studied the effect of process parameters to EG yield . Tai et al. (2012) reported that when 1.2% Ru/MC coupling was used with tungstic acid at a mass ratio of 2:1, maximum EG yield was achieved at 58.5% [2]. Wang and Zhang (2013) mentioned that by fine-tuning tungsten (W) content to transition metals, up to 60% EG yield could be attained [3]. Zheng et al. (2014) indicated that EG yield could be improved from 67% to 76% using more efficient tungstenic acid [4]. The work from Ji et al. (2008) showed that a maximum of 61% EG yield was possible when 2% Ni-30% W 2 C/AC-973 were used as catalyst [5] whereas using 5% Ni-15%W/SBA-15 and 10% Ni-30% W 2 C/AC, (both with mass ratio of Ni to W = 1:3), 76.1% and 73% EG yield could be achieved, respectively [6,7]. When Raney nickel to tungstic acid was used in a mass ratio of 2:3, maximum EG yield was achieved at 65.4% [8]. 61.0% EG yield was possible when 2% Ni-W2C/AC-973 was used [9]. Xiao et al. (2013) reported that a maximum of 31.6% EG yield was attained when CuCr (4) was coupled with 0.06 g Ca whereas with this system of catalysts at 10% cellulose loading, maximum EG yield of 7.6% was attained [10]. Using Jurusalem artichoke tuber, 14.1% EG yield was attained with 4% Ni-20% W 2 C/AC [11]. Using ammonium

•
It requires less resources (experiments, time, material, man hours, etc.) to get meaningful information; • The estimates of effects of each factor are more precise; • Interaction between factors can be estimated systematically; and • Experimental information in a larger region of the factor space could be obtained.
In addition, oil palm empty fruit bunches (EFB) fibres are regional-specific agricultural wastes, therefore none of the work published so far is able to predict EG yield when EFB is used as feedstock.
Since response surface methodology (RSM), specifically central composite design (CCD) has been widely used in biohydrogen [28] and biodiesel [29][30][31][32][33] production process optimization, it should therefore be used in this work to determine the effects of process parameters (binary catalyst ratio, temperature, pressure and tungstic acid to EFB ratio) of biomass hydrogenation to EG yield and eventually optimize EG yield.

Materials and Chemicals
All materials and chemicals were used as received unless stated otherwise. Empty fruit bunches (EFB) fibres were sourced from a local mill located in Selangor, Malaysia. The fibres were hydrothermally processed and hence free from oil and dirt, and were cut to a 1-inch length. Upon receiving, the fibres were further milled to 1-mm fibre length using a lab-scale grinder (IKA Labortechnik, Staufen, Germany) with a 1-mm sieve to get the desired fibre length of 1-mm and below. Raney Nickel was obtained from W.R. Grace and Company, Columbia, Maryland, USA. Analytical grade sodium hydroxide (NaOH) and tungstic acid were obtained from Sigma Aldrich (Subang Jaya, Selangor, Malaysia). Purified hydrogen was obtained from Mox-Linde, (Petaling Jaya, Malaysia).

Alkaline Pre-Treatment of Empty Fruit Bunches (EFB) Fibres
All EFB fibres used in this study had undergone alkaline pre-treatment in accordance to work reported by Zawawi et al. [34]. Alkaline pre-treatment of untreated EFB was carried out in a 1.5-L pressure vessel (PARR, Moline, Illinois, USA). 20 g of EFB and 180 g of 0.438 mol L −1 of NaOH solution at a ratio of one to nine parts of EFB to base solution were typically used. Alkaline pre-treatment of EFB was conducted at 423 K for 90 min at a stirring rate of 500 rpm. Upon cooling of the reactor, pre-treated EFB was obtained via filtration. The solid samples were washed repeatedly until neutralized pH was obtained. They were then dried in an oven until the moisture content was <10 wt. % and kept in a fridge until use. The characterization protocol of the pre-treated EFB was performed in accordance with the work reported by Zawawi et al. [34].

Catalytic Biomass Hydrogenation
Catalytic conversion of pre-treated EFB was carried out in a 1.5-L pressure vessel and 18.75-L pressure vessel (PARR, Moline, Illinois, USA). 20 g of pre-treated EFB and 180 g of deionized water at a ratio of one to nine parts of EFB to water were typically used in a 1.5-L pressure vessel, whereas 500 g of pre-treated EFB and 4500 g of deionized water were typically used in a 18.75-L pressure vessel. The mass ratio of Raney Nickel to tungstic acid varied from 0.80 to 1.20, temperature varied from 503 to 523 K, initial hydrogen pressure at ambient condition of 303 ± 2 K varied from 18 to 25 bar (g), whereas mass ratio of tungstic acid to EFB varied from 0.060 to 0.100. Residence time was fixed at 120 min and stirring rate was fixed at 500 rpm. After completion of catalytic biomass hydrogenation, the pressure vessel was depressurized, product gas was released as waste gas, product liquor was collected using vacuum filtration through 0.45 µm nylon filter, and proceed to EG yield quantification.

EG Yield Quantification
Polyhydric alcohol yield quantification was conducted using high-performance liquid chromatography (HPLC) (Agilent Technologies 1260 series) equipped with a refractive index detector (RID). Two (2) units of Shodex SC-1011 (8.0 mm I.D., 300 mm length) analytical columns connected in series, equipped with its guard column SC-LC (6.0 mm I.D., 50 mm length) (Showa Denko, Toyama-shi, Japan) were installed in a column holder with a column temperature set at 353 K. Deionized water was flowed at a flow rate of 0.4 mL min −1 and RID temperature was maintained at 313 K while 10 µL of product liquor was injected into the system.

Design of Experiment, Statistical Analysis and Optimization
In this study, process parameters binary catalyst ratio (Raney Nickel to tungstic acid, RN:TA), temperature, pressure and ratio of tungstic acid to EFB (TA:EFB) were represented with codes A, B, C and D. The relative contributions of process parameters to response EG yield were determined using RSM. A CCD was employed to investigate the optimal process parameters for EG yield, and 18 experiments with two replicated at the center point were carried out. A statistical software, Design-Expert ®® 8 (DX8 ®® ), developed by Stat-Ease, U.S.A. was used to generate randomized experimental conditions used for this work so that the effect of unexplained variability on the observed responses could be minimized. Table 1 summarizes the range and levels of each of the variables. The default CCD option for alpha set was to get a rotatable design with the axial (star) points set at 1.68179 coded unit from the center. This is a conventional choice for CCD and this should produce less extreme and hence a more practical value for our factors A-D [35].
Then RSM analysis was performed to obtain 3D surface graphs between variables (factors A, B, C, D) and response (EG yield). Analysis of variance (ANOVA) was used to statistically test the experimental results at a significance level of p = 0.05. It was also important to have a high p-value for Lack of Fit (LOF) before accepting the model. If p-value of LOF is low, i.e., ≤ 0.05, a check on residual plots is necessary to identify possible outliers and got them removed from the model in order to get an adequate model. Additionally, coefficient of determination (R 2 ) could be used to check the model adequacy.
Once model was ascertained, using the "Optimization" feature of DX8 ®® . Process parameters optimizing EG yield can be generated after the following is defined: • Determine goal of each factor and response to be either maximize, minimize or in range; • Set upper and lower limit for each factor and response; • Set importance weightage for each factor and response (1 being least important, 5 being paramount).
From the solutions generated by "Optimization", desirability closest to 1 were chosen to validate the model using the experiment.

Empty Fruit Bunches (EFB) Characterisation
The untreated and pre-treated EFB fibre characteristics were tabulated in Table 2. EFB characterization in this work wwas done in accordance to the method reported by Zawawi et al. [34]. Alkaline pre-treatment managed to give 54.31% delignification efficiency and comparatively increased cellulose content to 62.31% from its initial cellulose content before alkaline pre-treatment. These results further confirmed the EFB delignification model developed by Zawawi

Design of Experiments
Temperature (factor B) and hydrogen pressure (factor C) were chosen as parameters to study in this work because previously the work reported by Ji et al. and Luo et al. stated that they were the two main parameters affecting product yield and selectivity [9,13]. Moreover, hydrogen partial pressure would affect hydrogenation efficiency whereas temperature would affect retro-aldol reaction, and these two would eventually affect EG yield. The other two new parameters of mass ratio of binary catalyst (RN:TA) (factor A) as well mass ratio of catalyst to biomass (TA:EFB) (factor D) were chosen in this study due to their limited information reported. All 18 experiments run were completed and their respective EG yield was quantified. Factor A was 0.80-1.20, factor B was 18-25 bar (g), factor C was 230-250 • C (503-523 K) whereas factor D was 0.060-0.100. Table 3 summarizes the CCD experimental conditions and respective EG yield.
EG yield for all runs were in the range of 10.73 to 26.54 wt.%, with run no. 10 giving the lowest EG yield. When biomass hydrogenation was conducted at a temperature below 230 • C (503 K), more than 80% of biomass was found to be unreacted after the reaction ended. These results echoed the findings reported by Luo et al. that reaction temperatures above 503 K were required for stronger acidity that would favor cellulose hydrolysis, which was reflected by the disappearance of the partially hydrolysed products cellobiose and cellotriose [13]. Yu and Wu also reported that temperature played a vital role in hydrothermal conversion of cellulose, a temperature of 473-643 K and a pressure of 4-20 MPa were required to maintain water in its liquid state [43]. Ion product of water, H 3 O + and OH − , defined as the product of the concentration of acidic and basic forms of water, can be formed at this temperature range [43,44]. Massive ion product formation of water happened between 473 K and 643 K and achieved its maximum at 523 K [44,45]. This allowed water to produce H 3 O + protons at elevated temperatures, and promoted cellulose hydrolysis even without the presence of a catalyst [44].
At temperatures above 250 • C (523 K), sticky char formation was observed. The increase in temperature led to more water ion product formation, and higher acidity resulted in more C-C bond formation hence increasing the formation of sticky char. The remaining 11.08% lignin that was still present in the pre-treated EFB might have formed intramolecular condensation that led to sticky char formation. Through glycome profiling, Sebran et al. [46] discovered that 6.19 wt.% of lignin would still be present after harsh alkaline pre-treatment was performed on EFB, as alkaline pre-treatment could only remove the outer layer of lignin completely but the inner lignin structure would still remain intact.

Response Surface Analysis and Analysis of Variance (ANOVA)
The RSM model generated was first checked by examining the residuals plot and external studentized residuals plot, as shown in Figure 1A,B, respectively. Residuals showed the difference between actual and predicted response values [35]. Assuming the model holds, external studentized residual was used to check whether a run was consistent with other runs, whereas the model coefficients were calculated based on all design points except one and a prediction of the response at this point was made, where any value greater than the calculated standard error, usually between absolute value of 3 and 4, meant that this point should be examined as a possible outlier [35].
Referring to Figure 1B, run no. 10 appeared to be an outlier as it estimated standard errors of −4.752, which was larger than the calculated standard error of ±3.748. Figure 1A further confirmed that run no. 10 was an outlier due to its large difference between actual and predicted values of −6.71. This explained why an EG yield of run no. 10 showed abnormal results of lower-than-normal EG yield. Therefore, run no. 10 was subsequently ignored and new plots were generated in Figure 2 showing no outliers in these plots. Table 4 shows the fit summary of the response surface model upon removing the outlier. The linear model and quadratic model were suggested by DX8 ®® due to their model p-value that was below 0.05. Although the linear model had a lower p-value of 0.0019, its lower p-value of 0.0241 for lack of fit (LOF) and lower adjusted R 2 value of 0.6455 were however not desirable for a good model. Hence, the quadratic model with a p-value of 0.0227 but higher adjusted R 2 value of 0.9845 and higher p-value of 0.0697 for LOF, would be desirable for our application. Selection criteria is the model with the highest order of polynomial where the additional terms are significant, and the model should not be aliased [35].
In order to obtain a more significant model and improve LOF (by having higher p-value of LOF), ANOVA was reduced by removing one insignificant factor, i.e., A 2 (the second-order term confounded with main factor A). ANOVA of the reduced quadratic model and its fit statistics are summarized in Table 5.  The model had an F-value of 112.13, indicating the model was significant. There was only a 0.12% chance that an F-value this large could occur due to noise. The Lack of Fit (LOF) of this model had an F-value of 43.88, implying that the LOF was not significant relative to the pure error. There was a 10.61% chance that a LOF with the F-value this large could occur due to noise. An insignificant LOF was good as we wanted the model to fit [35]. Insignificant LOF indicates that our model does not have unexplainable differences between replicates [35]. The reduced quadratic model had an adjusted R 2 value of 0.9890, and this adjusted R 2 calculated the percentage of response variation contributed by significant factors, i.e., factors with a p-value < 0.05. Adequate Precision measures the signal to noise ratio and a ratio greater than 4 is desirable for a good prediction model [35]. This reduced quadratic model had an adequate precision of 38.885, and therefore it could be used to navigate the design space and also to optimize EG yield by varying the four factors studied [35].  Table 4 shows the fit summary of the response surface model upon removing the outlier. The linear model and quadratic model were suggested by DX8 ®® due to their model p-value that was below 0.05. Although the linear model had a lower p-value of 0.0019, its lower p-value of 0.0241 for lack of fit (LOF) and lower adjusted R 2 value of 0.6455 were however not desirable for a good model. Hence, the quadratic model with a p-value of 0.0227 but higher adjusted R 2 value of 0.9845 and higher p-value of 0.0697 for LOF, would be desirable for our application. Selection criteria is the model with the highest order of polynomial where the additional terms are significant, and the model should not be aliased [35].  The ANOVA summary showed that EG yield was highly dependent on all four factors studied, i.e., A, B, C, D and second order terms confounded with main factors like AC, AD, BC, BD, CD, B 2 , C 2 , D 2 were also significant model terms. Figure 3 shows the factor interaction contour between A and C Processes 2020, 8, 588 9 of 13 as well as C and D, with respect to EG yield. Operating factors that led to the red zone would lead to the desired outcome of higher EG yield of above 25 wt.%.
The ANOVA summary showed that EG yield was highly dependent on all four factors studied, i.e., A, B, C, D and second order terms confounded with main factors like AC, AD, BC, BD, CD, B 2 , C 2 , D 2 were also significant model terms. Figure 3 shows the factor interaction contour between A and C as well as C and D, with respect to EG yield. Operating factors that led to the red zone would lead to the desired outcome of higher EG yield of above 25 wt.%.

EG Yield Optimization and Verification
By using this model and employing "Optimization" feature of DX8 ®® , we set constraints for each factor which included their corresponding goal, limit and importance, in order to maximize EG yield, as shown in Table 6. By applying the desirability function method, the optimal process parameters that gave a desirability of 1 suggested optimal EG yield to be at 26.89 wt.% with the following conditions:
Pressure vessel with different sizes (1.5-L and 18.75-L) were used to verify the suggested optimal process parameters and Table 7 shows the results of these verification runs. Under the optimal conditions, the predicted and actual EG yield were very close and within the standard deviation of the data. The suggested optimum EG yield from this quadratic model is >25 wt.% and our RSM model was successfully ascertained and reproduced in the laboratory. This quadratic model could be used to predict EG yield for future work.

Conclusions
Using DX8 ®® , a quadratic model was developed which identifies factors like binary catalyst ratio (A), pressure (C) and catalyst to feedstock ratio (D), and their confounding factors as significant. Temperature (B) was comparatively less significant due to the narrow temperature range that we studied. As suggested by DX8 ®® , one of the optimum parameters to obtain an optimized EG yield of 26.89 wt.% would be having a binary catalyst ratio (Raney Nickel to tungstic acid) of 0.96, a temperature of 240 • C (513 K), a pressure of 24.8 bar (g) and catalyst (tungstic acid) to EFB ratio of 0.10.
Having a high coefficient of determination (R 2 ) value of more than 0.98, the quadratic model developed in this study could be used to give a precise EG yield prediction even when the four parameters are being modified in the future. Although many research groups had reported their findings on glucose, cellulose or biomass hydrogenation for EG production, there wasn't any literature reporting effects of factors as well as effects of confounding factors affecting EG yield. This study has been able to fill in this missing gap nicely.