Multiobjective Optimization for the Greener Synthesis of Chloromethyl Ethylene Carbonate by CO 2 and Epichlorohydrin via Response Surface Methodology

: In this paper, a statistical analysis with response surface methodology (RSM) has been used to investigate and optimize process variables for the greener synthesis of chloromethyl ethylene carbonate (CMEC) by carbon dioxide (CO 2 ) and epichlorohydrin (ECH). Using the design expert software, a quadratic model was developed to study the interactions e ﬀ ect between four independent variables and the reaction responses. The adequacy of the model was validated by correlation between the experimental and predicted values of the responses using an analysis of variance (ANOVA) method. The proposed Box-Behnken design (BBD) method suggested 29 runs for data acquisition and modelling the response surface. The optimum reaction conditions of 353 K, 11 bar CO 2 pressure, and 12 h using fresh 12% ( w / w ) Zr / ZIF-8 catalyst loading produced 93% conversion of ECH and 68% yield of CMEC. It was concluded that the predicted and experimental values are in excellent agreement with ± 1.55% and ± 1.54% relative errors from experimental results for both the conversion of ECH and CMEC yield, respectively. Therefore, statistical modelling using RSM can be used as a reliable prediction technique for system optimization for greener synthesis of chloromethyl ethylene carbonate via CO 2 utilization.


Introduction
Carbon dioxide (CO 2 ) chemistry has earned enormous interest in recent years due to its abundance and inexpensive nature. It is a nontoxic, non-flammable, easily available, and typical renewable C1 source of organic synthesis [1]. CO 2 is an important "greenhouse" gas that has drawn greater attention in line with the need for the development of green engineering and sustainable society. In this regard, the development of environmentally benign and efficient synthetic of chemical utilization of CO 2 has been a subject of immense research in academia as evidenced by the rising number of publications in all areas of CO 2 management [2]. Although CO 2 fixation is unlikely to consume large quantities of CO 2 in the atmosphere, this measure can be regarded as a significant strategy for the development of sustainable and safe processes [3]. With the intriguing applications of organic carbonates, the use of CO 2 as a raw material to synthesize cyclic organic carbonates has gained extensive attention in chemical industries [4].
Organic carbonates are versatile compounds used as raw materials for many industrial applications including raw materials for polycarbonates and polyurethane synthesis [5], green solvents [6], techniques, which are full three-level factorial designs: Box-Behnken designs, central composite designs, and Doehlert designs.
A multivariate optimization technique is a statistical tool for analyzing complex non-linear processes. This is especially useful when interactions are not known or optimal process parameters are to be determined in order to make a process more robust [43]. It is cost-effective as fewer experimental trials are required, high computational efficiency [44], and it requires very little or no human experience to obtain an accurate and satisfactory results [45]. Therefore, the systematic application of RSM optimization for the catalytic conversion of epichlorohydrin (ECH) and carbon dioxide (CO 2 ) to chloromethyl ethylene carbonate (CMEC) can be regarded as an innovative way of CO 2 utilization.

Catalysts Preparation
Zirconium-doped ZIF-8 (Zr/ZIF-8) was synthesized according to a method, which was previously described elsewhere [46,47]. Briefly, 8 mmol of zinc nitrate hexahydrate (Zn(NO 3 ) 2 ·6H 2 O, purity 99.99%) and zirconium (IV) oxynitrate hydrate (ZrO(NO 3 ) 2 ·6H 2 O, purity 99.99%) solutions in a stoichiometric ratio of Zn: Zr = 9:1 were dissolved in 6.2 mmol of methanol. A separate solution of 14.2 mmol of 2-methylimidazole and 600 mml of methanol was prepared in another flask which was added by dropwise addition to the Zr-Zn based solution. The mixture, conducted in an ambient temperature under nitrogen gas flow, was vigorously stirred for 6 hrs. The crystals were collected and separated by centrifugation at 300 rpm for 30 min. The solution was washed thoroughly with methanol three times and then dried at room temperature. The crystals were left to dry overnight at 373 K. The greyish-white powders of Zr-ZIF-8 sample were further washed with DMF for 24 h in order to remove any excess of an unreacted organic linker. The solution was then heated at a temperature of 373 K in order to activate it. The sample was allowed to cool down to room temperature naturally before been capped in a vial and refrigerated, which was ready for use in catalytic reactions.

Proposed Reaction Mechanism and Reaction Pathways
On the basis of our experimental results and theoretical understanding, we proposed a plausible reaction mechanism for the coupling reaction of ECH and CO 2 . Figure 1 shows the reaction mechanism was initiated by coordination of ECH with Lewis acid site Zn 2+ to form the adduct of zinc-epoxide complex, then nucleophilic interaction on the electrophilic carbon of CO 2 (step 1). At the same time, the acidic sites (unsaturated coordinative Zn or structural defects) of Zr/ZIF-8 interact with the oxygen atom of an epoxide (step 2). The activated CO 2 attacks the less sterically hindered carbon atom of epoxide, which results in the epoxide ring-opening (step 3). Finally, the ring-closure step takes place between the O−anion and carbon atom in the intermediates to produce CMEC (step 4). Figure 2 shows the reaction pathways 1, 2, and 3 with some by-products. The decline in selectivity and CMEC yield was expected because the gas chromatography mass spectroscopy (GC-MS) analysis of the samples shows that 17.3% of 3-chloropropane 1,2-diol and 14.1% 2,5-bis (chloromethyl)-1,4-dioxane (by-products) have been formed at 353 K. Similar by-products and results have been reported by Zhou et al. [48]. This may explain in part why a drop in selectivity and yield of CMEC was recorded.

One-Factor-at-a-Time (OFAT) Analysis
OFAT analysis was developed to determine the preliminary effective range of the selected parameters for statistical analysis. The effect of four-single factors (temperature, pressure, reaction time, and catalyst loading) were evaluated for the synthesis of chloromethyl ethylene carbonate. The OFAT analysis investigated all the four parameters in the following range: Reaction

Experimental Design
Based on the OFAT results, a three-level, four-factor (3 4 ) factorial design with 29 runs of experiments were suggested for this study in order to determine the responses (conversion and yield). In this design, all the four factors were varied simultaneously over a set of experimental runs. To avoid bias, the suggested set of experiments were carried out randomly and the four factors, temperature, pressure, catalyst loading, and reaction time have been labelled as x 1 , x 2 , x 3 , and x 4 , respectively as shown in Table 1. The variables and their coded and uncoded values are presented with each level and range as given below in Table 1 (i.e., −1, 0, 1). The total number of experiments (N) is given by Equation (1): where, k is the number of independent variables, Cp is the replicate number of the center point.

Statistical Analysis
The empirical mathematical model showing the effect of the independent variables x 1 , x 2 , x 3 , and x 4 on the predicted response Y was investigated using the second order polynomial regression equation with backward elimination.
A quadratic equation derived using RSM for the model is shown using Equation (2): where Y is the predicted response, x i and x j are the independent variables in coded levels (I j), b i , b ii , and b ij are the coefficients for linear, quadratic, and interaction effects, respectively, b 0 is the model coefficient constant, n is the number of factors, and ε is the model random error [49]. The adequacy of the predicted models was validated by a number of statistical tools such as correlation coefficient (R 2 ), adjusted coefficient of determination (R 2 adj ), and the predicted coefficient of determination (R 2 pred ). The statistical significance of the predicted model was analyzed by (ANOVA) using a regression coefficient by conducting the Fisher's F-test at 95% confidence level [50]. Design Expert 11 software (Stat-Ease Inc., Minneapolis, MN, USA) was used for the design of experiment, regression, and graphical analysis. Statistical significance of the results have been presented by p < 0.05 and mean ± SE. The fit quality of the polynomial equation has been proved by R 2 .

Experimental Procedures
In a typical cycloaddition reaction, chloromethyl ethylene carbonate (CMEC) was synthesized from epichlorohydrin (ECH) and carbon dioxide (CO 2 ) in a solvent free and co-catalyst free conditions. A 25 mL stainless steel autoclave reactor equipped with a stirrer, thermocouple, heating mantle, and controller was initially charged with the required amount of limiting reactant ECH and a known amount of Zr/ZIF catalyst. The reactor was then heated to a specific temperature and continuously stirred. When the desired reaction temperature was reached, a known amount of liquid CO 2 was injected to the reactor via SCF pump at an assumed t = 0. The reaction mixture was left stirring and monitored for a set period of time.
After the reaction was completed, the reactor was cooled down to room temperature using an ice bath, depressurized, and then the reaction mixture was filtered. The catalyst was separated, washed with acetone, and dried in a vacuum oven. The product obtained from the filtered reaction mixture was then analyzed using gas chromatography (GC).

Analysis of Variance (ANOVA)
An analysis of variance (ANOVA) was performed using the Design Expert software in order to investigate the fitness and significance of the model for each regression coefficient. The empirical analysis of RSM model used to correlate the interactive relationship between the controlling factors (x 1 , x 2 , x 3 , and x 4 ) and the predicted response Y (conversion of ECH and yield of CMEC) are shown in Table 2 above. The results of the experimental trials at various process conditions show the range of the responses from 42% to 93% of ECH conversion and 16% to 68% of CMEC yield. This trend is consistent with the results published by Onyenkeadi and colleagues. The predicted values sufficiently correlate with the observed values and fit the RSM model design for this study. The best fitting model was established by a regression analysis using Design Expert software. Fitting of the data to various models (linear, two factors interactions (2FI), quadratic, and cubic polynomials) and their following analysis of variance (ANOVA).

Development of Regression Model
In this study, the purpose of using the RSM was to generate a statistical model that demonstrated mutual interaction between the responses and the effective variables. Through the experimental matrix generated in a randomized run of experiments, the obtained responses are given using second order polynomial regression equation with backward elimination as shown below. The equations show the empirical relationship between the conversion of ECH and the yield of CMEC and the experimental factors in coded forms. Y 1 and Y 2 are the response variables: ECH conversion and CMEC yield. The independent variables are x 1, x 2 , x 3 , and x 4 which are reaction temperature, pressure, catalyst loading, and reaction time, respectively. The results of interaction effects between the independent variables were deduced as follows: Temperature-pressure; x 1 x 2 , temperature-catalyst loading; x 1 x 3 , temperature-reaction time; x 1 x 4 , pressure-catalyst loading; x 2 x 3 , pressure-time; x 2 x 4 and catalyst loading-reaction time; x 3 x 4 . Finally, the excess of each independent variable was represented as follows: Temperature-temperature; x 1 2 , pressure-pressure; x 2 2 , catalyst loading-catalyst loading; x 3 2 and reaction time -reaction time; x 4 2 .

Statistical Analysis of Regression Model
The response model calculated for this study has demonstrated a high degree of accuracy with an R 2 of 0.9973 and an R 2 adj of 0.9954 at a confidence level of 95%. This agrees well with the result of [51] where the determination coefficient values, R 2 and R 2 adj , for the reliability of the model fitting, were calculated to be 0.9932 and 0.9658, respectively. Mäkelä et al. [52], also suggested that a good model fit should yield an R 2 of at least 0.8. Furthermore, the values of R 2 and R 2 adj are close to 1.0. This demonstrates that a mutual correlation exists between the experimental and the predicted values. Therefore, the statistical significance of the second-order polynomial equation for this design shows that the regression model is statistically significant (p < 0.0001) and the lack of fit test is non-significant (p > 0.05) relative to the pure error.
The following assumptions have been used to conclude the statistical adequacy checking of the model based on the ANOVA results. The first assumption is the similarity between the predicted and actual data of the two models as shown in Figure 3. This demonstrated that the variations between Figure 3a,b are statistically non-significant (NS) and the predicted model can be said to show a high level of accuracy and adequacy. Another assumption is the normality of the residuals. The plot of residuals has been investigated using normal plot where most of the points approximately form a straight line as shown in Figure 4a,b. This shows that residuals for both ECH conversion and CMEC yield are in normal distribution. This assumption is consistent with the report of Mäkelä et al. [52]. Thirdly, the randomization of the residuals have also been assessed using a plot between the residuals versus predicted responses. The random distribution in Figure 5 shows lack of clear structure with a normal distribution at zero mean and variance [53]. It can be observed in Figure 5a,b that points above and below the diagonal line show areas of over or under prediction with no definite structure.

Model Fitting and Adequacy Checking
In order to verify the model for fitting and adequacy test at 95% confidence level, it was necessary to apply analysis of variance (ANOVA). As shown in Tables 3 and 4, the ANOVA results indicated a good model fit with the model. F-value is 0.44 (Table 3) and the probability > F of less than 0.0001 implied that this model was significant. The lack of fit test (non-significant: p > 0.05) was also considered a good statistical indicator for the model adequacy checking as it relates the residual error to the pure error from the replica design point [54]. As indicated in ANOVA Tables 3 and 4, the conversion of ECH and CMEC yield was significantly (p < 0.05) influenced by the interactive and quadratic effects of all the independent variables.

Response Surface Plots Analysis
After the regression models had been built and model adequacy checking was tested, 3D response surface plots and their corresponding 2D contour plots were drawn for a model equation. Different shapes of the contour plots indicate different levels of interaction between two variables. For example, an oval plot represents significant interactions between the two selected variables while a circular plot means otherwise [54]. According to Rabiee et al. [55], 3D response surface promotes understanding of system behavior. It is also significant in recognizing the characters of response surface [56].

Effect of One Factor at a Time Experiments on Responses (OFAT)
The effects of individual reaction variables (temperature, pressure, time, and catalyst loading) and their interactions on reaction responses (conversion and yield) have been investigated using the 3D-surface and 2D-contour plots generated from the predicted quadratic model as evidenced in Figures 6-9. The experiments have been carried out by varying one reaction parameter at a time while keeping other parameters constant at the following reaction conditions: Reaction temperature 353 K, CO 2 pressure 11 bar, reaction time 12 h, catalyst loading 12% (w/w).

Response Surface Plots Analysis
After the regression models had been built and model adequacy checking was tested, 3D response surface plots and their corresponding 2D contour plots were drawn for a model equation. Different shapes of the contour plots indicate different levels of interaction between two variables. For example, an oval plot represents significant interactions between the two selected variables while a circular plot means otherwise [54]. According to Rabiee et al. [55], 3D response surface promotes understanding of system behavior. It is also significant in recognizing the characters of response surface [56].

Effect of One Factor at a Time Experiments on Responses (OFAT)
The effects of individual reaction variables (temperature, pressure, time, and catalyst loading) and their interactions on reaction responses (conversion and yield) have been investigated using the 3D-surface and 2D-contour plots generated from the predicted quadratic model as evidenced in Figures 6-9. The experiments have been carried out by varying one reaction parameter at a time while keeping other parameters constant at the following reaction conditions: Reaction temperature 353 K, CO2 pressure 11 bar, reaction time 12 h, catalyst loading 12% (w/w).

Response Surface Plots Analysis
After the regression models had been built and model adequacy checking was tested, 3D response surface plots and their corresponding 2D contour plots were drawn for a model equation. Different shapes of the contour plots indicate different levels of interaction between two variables. For example, an oval plot represents significant interactions between the two selected variables while a circular plot means otherwise [54]. According to Rabiee et al. [55], 3D response surface promotes understanding of system behavior. It is also significant in recognizing the characters of response surface [56].

Effect of One Factor at a Time Experiments on Responses (OFAT)
The effects of individual reaction variables (temperature, pressure, time, and catalyst loading) and their interactions on reaction responses (conversion and yield) have been investigated using the 3D-surface and 2D-contour plots generated from the predicted quadratic model as evidenced in Figures 6-9. The experiments have been carried out by varying one reaction parameter at a time while keeping other parameters constant at the following reaction conditions: Reaction temperature 353 K, CO2 pressure 11 bar, reaction time 12 h, catalyst loading 12% (w/w).

Effect of Reaction Temperature
To a significant extent, it is largely agreed that a directly proportional relationship exists between temperature and CMEC yield as shown in the results of ANOVA in Table 4. The influence of reaction temperature on CMEC yield has been investigated by varying temperature over the range of 323 K to 373 K. As evidenced in Figure 6, CMEC yield increased steadily from 40% to 68% as temperature increased from 323 K to 353 K. However, a gradual decrease in CMEC yield was observed at higher temperature values beyond 353 K. This may be due to the formation of diols and dimers of epichlorohydrin above optimum temperature [57]. Saada et al. and co-workers explained that higher reaction temperatures caused a shift in the equilibrium to the reactant side and resulted in a reduced DMC yield. The same temperature effect was also reported by Kilic et al. [58]; they have observed that as they increased the reaction temperature from 348 K to 373 K (while keeping other variables constant), there was a corresponding increase in ECHC yield from 65.8% to 97.0%. However, further increase in temperature beyond 373 K, caused a slight decrease both in the ECHC yield and catalyst selectivity.

Effect of Reaction Temperature
To a significant extent, it is largely agreed that a directly proportional relationship exists between temperature and CMEC yield as shown in the results of ANOVA in Table 4. The influence of reaction temperature on CMEC yield has been investigated by varying temperature over the range of 323 K to 373 K. As evidenced in Figure 6, CMEC yield increased steadily from 40% to 68% as temperature increased from 323 K to 353 K. However, a gradual decrease in CMEC yield was observed at higher temperature values beyond 353 K. This may be due to the formation of diols and dimers of epichlorohydrin above optimum temperature [57]. Saada et al. and co-workers explained that higher reaction temperatures caused a shift in the equilibrium to the reactant side and resulted in a reduced DMC yield. The same temperature effect was also reported by Kilic et al. [58]; they have observed that as they increased the reaction temperature from 348 K to 373 K (while keeping other variables constant), there was a corresponding increase in ECHC yield from 65.8% to 97.0%. However, further increase in temperature beyond 373 K, caused a slight decrease both in the ECHC yield and catalyst selectivity.

Effect of Reaction Temperature
To a significant extent, it is largely agreed that a directly proportional relationship exists between temperature and CMEC yield as shown in the results of ANOVA in Table 4. The influence of reaction temperature on CMEC yield has been investigated by varying temperature over the range of 323 K to 373 K. As evidenced in Figure 6, CMEC yield increased steadily from 40% to 68% as temperature increased from 323 K to 353 K. However, a gradual decrease in CMEC yield was observed at higher temperature values beyond 353 K. This may be due to the formation of diols and dimers of epichlorohydrin above optimum temperature [57]. Saada et al. and co-workers explained that higher reaction temperatures caused a shift in the equilibrium to the reactant side and resulted in a reduced DMC yield. The same temperature effect was also reported by Kilic et al. [58]; they have observed that as they increased the reaction temperature from 348 K to 373 K (while keeping other variables constant), there was a corresponding increase in ECHC yield from 65.8% to 97.0%. However, further increase in temperature beyond 373 K, caused a slight decrease both in the ECHC yield and catalyst selectivity.

Effect of Carbon Dioxide (CO 2 ) Pressure
ANOVA Table 4 demonstrates the dependence of CO 2 pressure on CMEC yield, since CO 2 acts both as reactant and reaction medium simultaneously [59]. As indicated in Figure 7, when CO 2 pressure was increased from 8 to 11 bar, the CMEC yield also increased from 50% to 68%. Conversely, with the CO 2 pressure of 11.5 bar, a 59% CMEC yield was recorded indicating a declining effect. Zhong et al. [60] demonstrated the effect of variation in CO 2 pressure on organic carbonates. They have enhanced more propylene carbonate (PC) yield when CO 2 pressure was increased from 1 MPa to 3 MPa. However, when CO 2 pressure was further increased to 4 MPa, they observed that the concentration of propylene oxide (PO) in gas phase had decreased as a result of dilution by CO 2 and consequently resulted in a reduced PC yield. It is therefore concluded that the optimum CO 2 pressure based on OFAT analysis for this set of experiments was 11 bar of CO 2 pressure.

Effect of Reaction Time
Reaction time is one of the crucial factors in a catalytic reaction. Figure 8 shows a direct proportionality effect between reaction time and the CMEC yield; the yield increased gradually as reaction time increased until it reached 68% in 12 h. Further increase in reaction time beyond 12 h resulted in a continuous decline in CMEC yield as shown in Figure 8. This could be as a result of formation of polymerized CMEC caused by prolonged reaction time [61]. A similar phenomenon was also reported by Onyenkeadi and co-workers, where increase in reaction time from 8 to 16 h was directly proportional to butylene carbonate (BC) yield. However, prolonged reaction time beyond this time resulted in decrease in BC yield.

Effect of Catalyst Loading
The effect of catalyst loading on CMEC yield was investigated by varying Zr/ZIF-8 loading from 5% to 15% (w/w). As shown in Figure 9, it can be observed that as catalyst loading was increased, CMEC yield also increased proportionally from 42% reaching a maximum of 68% at 12% (w/w) catalyst loading. It was then decreased progressively when the amount of catalyst was further increased to 13% (w/w), indicating that optimum catalyst loading had been exceeded. It would be expected that the number of active sites available for the reaction of ECH and CO 2 would increase as catalyst loading increased [62]. However, Han et al. [63] argued that an excessive increase in catalyst loading tends to provoke formation of undesirable side-products (in their experiment, a by-product of diglyceride (GDL) or triglyceride (GTL) was formed), thereby causing a drop in glycerol monolaurate (GML) selectivity as they increased the amount of catalyst beyond 2% (w/w). Similarly, in the present work, increase in the amount of catalyst loading beyond the optimum level was unfavorable to the reactive system resulting in a reduced CMEC yield. Therefore, the optimum catalyst loading for this reactive system is 12% at a reaction temperature of 353 K for 12 h at 11 bar of CO 2 pressure.

Interactive Effect of Process Variables on Responses
The interaction effect of each pair of reaction variables have been investigated using ANOVA results, 3D surface, and 2D contour plots. The interaction effect of some process variables on ECH conversion and CMEC yield produce different effects at different levels of other variables. Therefore, 3D plots have played a crucial role in making accurate predictions about process optimization [64]. From ANOVA Tables 3 and 4, it can be observed that all the four reaction parameters are deemed significant and can affect the process response tremendously at different levels of interaction. Hence, the interactive effect of process variables has a direct influence on the system optimization. The interaction effect between a pair of variables would be negligible if the contour plot of the response surface is circular. Conversely, the interactions effect would be significant if the contour plot is elliptical [65]. Therefore, instead of studying single variable (as in conventional method) the interactions were investigated which is significant for a comprehensive optimization study.

Interactive Effect of Temperature and Pressure
As depicted in Figure 10 and the ANOVA Tables 3 and 4, the interaction effect of reaction temperature and CO 2 pressure has played significant roles in both ECH conversion and CMEC yield (while keeping reaction time and catalyst loading at their optimum: 12 h and 12% (w/w), respectively). At lower reaction temperature (e.g., at 323 K), increase in the CO 2 pressure from 4 to 16 bar increases the CMEC yield from 47% to 68%. However, higher reaction temperature beyond 353 K showed a negative effect on CMEC yield ( Figure 10a); this could possibly be as a result of formation of by-products at elevated temperature as indicated in the reaction mechanism ( Figure 1). Furthermore, at a different level of interaction between temperature and pressure (e.g., from 358 K to 373 K and 13-16 bar), a notable effect was also recorded where there was a gradual decline in the CMEC yield indicating optimum condition had been exceeded. This shows that variation in reaction temperature had a negative effect on both responses at higher values. Therefore, the temperature-pressure relationship has significant effect on process responses. Similarly, the elliptical shape of the 2D contour plot in Figure 10b exemplifies a mutual interactive effect of the reaction variables on responses.

Interactive Effect of Temperature and Pressure
As depicted in Figure 10 and the ANOVA Tables 3 and 4, the interaction effect of reaction temperature and CO2 pressure has played significant roles in both ECH conversion and CMEC yield (while keeping reaction time and catalyst loading at their optimum: 12 h and 12% (w/w), respectively). At lower reaction temperature (e.g., at 323 K), increase in the CO2 pressure from 4 to 16 bar increases the CMEC yield from 47% to 68%. However, higher reaction temperature beyond 353 K showed a negative effect on CMEC yield (Figure 10a); this could possibly be as a result of formation of byproducts at elevated temperature as indicated in the reaction mechanism ( Figure 1). Furthermore, at a different level of interaction between temperature and pressure (e.g., from 358 K to 373 K and 13-16 bar), a notable effect was also recorded where there was a gradual decline in the CMEC yield indicating optimum condition had been exceeded. This shows that variation in reaction temperature had a negative effect on both responses at higher values. Therefore, the temperature-pressure relationship has significant effect on process responses. Similarly, the elliptical shape of the 2D contour plot in Figure  10b Figure 11 illustrates the interaction effect of reaction time and temperature on CMEC yield (while keeping other two variables at their optimum: Catalyst loading: 12% (w/w), CO 2 pressure: 11 bar). The surface plot suggested that the CMEC yield was highest (68%) at a reaction time of 12 h and temperature of 353 K, indicating that an increase in the reaction temperature from 313 K to 353 K favors ECH conversion and consequently enhances CMEC yield as shown in Figure 11a. However, increase in reaction temperature beyond 353 K at 12 h of reaction time was unfavorable to the reactive system causing a marginal drop in CMEC yield from 68% to 65%. Onyenkeadi et al. and co-workers reported that formation of oligomers and isomers are possible at extended reaction time at higher temperature. Product quality and stability may also be affected due to chemical degradation or losses by thermal decomposition at higher reaction temperature [66]. Response surface and contour plots of Figure 11 clearly show that CMEC yield had a linear effect with increasing reaction temperature until the optimum condition was achieved. This phenomenon agrees with the Arrhenius law [67]; higher temperature results in a higher conversion rate and consequently leading to higher CMEC yield. It can be concluded from the ANOVA Table 4 that the reaction temperature was found to be a highly influencing parameter on both the conversion of ECH and CMEC yield as evident from low p-value (<0.0001).

Interactive Effect of Temperature and Time
Energies 2020, 13, x FOR PEER REVIEW 16 of 27 Figure 11 illustrates the interaction effect of reaction time and temperature on CMEC yield (while keeping other two variables at their optimum: Catalyst loading: 12% (w/w), CO2 pressure: 11 bar). The surface plot suggested that the CMEC yield was highest (68%) at a reaction time of 12 h and temperature of 353 K, indicating that an increase in the reaction temperature from 313 K to 353 K favors ECH conversion and consequently enhances CMEC yield as shown in Figure 11a. However, increase in reaction temperature beyond 353 K at 12 h of reaction time was unfavorable to the reactive system causing a marginal drop in CMEC yield from 68% to 65%. Onyenkeadi et al. and co-workers reported that formation of oligomers and isomers are possible at extended reaction time at higher temperature. Product quality and stability may also be affected due to chemical degradation or losses by thermal decomposition at higher reaction temperature [66]. Response surface and contour plots of Figure 11 clearly show that CMEC yield had a linear effect with increasing reaction temperature until the optimum condition was achieved. This phenomenon agrees with the Arrhenius law [67]; higher temperature results in a higher conversion rate and consequently leading to higher CMEC yield. It can be concluded from the ANOVA Table 4 that the reaction temperature was found to be a highly influencing parameter on both the conversion of ECH and CMEC yield as evident from low p-value (<0.0001).

Interactive Effect of Temperature and Catalyst Loading
The overall CMEC yield has been significantly influenced by the interaction between the catalyst loading and reaction temperature while CO 2 pressure and time have been kept at optimum values of 11 bar and 12 h, respectively. For example, Figure 11 shows that at lower catalyst loading of 5% (w/w), only 34% of CMEC yield was recorded as a result of low ECH conversion at low catalyst loading. The CMEC yield increased steadily up to 68% as reaction temperature increased at moderate levels of catalyst loading from 333 K to 353 K. This phenomenon could be attributed to the increase in the catalyst surface area, which provides more contact area between the limiting reactant ECH and the active sites of the catalyst. Higher catalyst loading gives higher ECH conversion resulting in higher CMEC yield-an effect which is more pronounced at higher temperatures. However, at higher temperature above 353 K, a marginal decrease in CMEC yield was observed, which may be due to catalyst deactivation at very high temperature [68]. The contour plot in Figure 12b with elliptical shape demonstrated the significant and combined effect of the catalyst loading and reaction temperature. The result has also supported lower p-value (0.0005) of the interaction x 1 x 3 term. As shown in Figure 12a, at any designated value of reaction temperature from 333 K to 353 K, the CMEC increased proportionally with catalyst loading. This observation was also supported by low p-value (<0.0001).

Interactive Effect of Temperature and Catalyst Loading
The overall CMEC yield has been significantly influenced by the interaction between the catalyst loading and reaction temperature while CO2 pressure and time have been kept at optimum values of 11 bar and 12 h, respectively. For example, Figure 11 shows that at lower catalyst loading of 5% (w/w), only 34% of CMEC yield was recorded as a result of low ECH conversion at low catalyst loading. The CMEC yield increased steadily up to 68% as reaction temperature increased at moderate levels of catalyst loading from 333 K to 353 K. This phenomenon could be attributed to the increase in the catalyst surface area, which provides more contact area between the limiting reactant ECH and the active sites of the catalyst. Higher catalyst loading gives higher ECH conversion resulting in higher CMEC yield-an effect which is more pronounced at higher temperatures. However, at higher temperature above 353 K, a marginal decrease in CMEC yield was observed, which may be due to catalyst deactivation at very high temperature [68]. The contour plot in Figure 12b with elliptical shape demonstrated the significant and combined effect of the catalyst loading and reaction temperature. The result has also supported lower p-value (0.0005) of the interaction x1x3 term. As shown in Figure 12a, at any designated value of reaction temperature from 333 K to 353 K, the CMEC increased proportionally with catalyst loading. This observation was also supported by low p-value (<0.0001).

Interactive Effect of Time and Pressure
Similar to the previous observation of the interaction effect of temperature and pressure, Figure 13 demonstrates the interaction effect of CO 2 pressure and time on CMEC yield while maintaining reaction temperature and catalyst loading at 353 K and 12% (w/w), respectively. For example, at a shorter reaction time of 4 h, there was a negligible effect of CO 2 pressure in the CMEC yield. Figure 13a shows that optimum reaction time of 12 h was observed at a CO 2 pressure of about 12 bar with a 68% of CMEC yield. It has been observed in Figure 13 that the CMEC yield reached a maximum at a reaction time of 12 h, thereafter, it was stable. A further increase in reaction time beyond this value caused a sharp drop in CMEC yield as indicated in surface plot of Figure 13b.

Interactive Effect of Time and Pressure
Similar to the previous observation of the interaction effect of temperature and pressure, Figure 13 demonstrates the interaction effect of CO2 pressure and time on CMEC yield while maintaining reaction temperature and catalyst loading at 353 K and 12% (w/w), respectively. For example, at a shorter reaction time of 4 h, there was a negligible effect of CO2 pressure in the CMEC yield. Figure 13a shows that optimum reaction time of 12 h was observed at a CO2 pressure of about 12 bar with a 68% of CMEC yield. It has been observed in Figure 13 that the CMEC yield reached a maximum at a reaction time of 12 h, thereafter, it was stable. A further increase in reaction time beyond this value caused a sharp drop in CMEC yield as indicated in surface plot of Figure 13b.

Interactive Effect of Time and Catalyst Loading
The contour and 3D surface plots in Figure 14 show the interaction effect between the reaction time and the catalyst loading at a constant temperature of 353 K and CO2 pressure of 12 bar. The

Interactive Effect of Time and Catalyst Loading
The contour and 3D surface plots in Figure 14 show the interaction effect between the reaction time and the catalyst loading at a constant temperature of 353 K and CO 2 pressure of 12 bar. The contour plots show less curvature up to 7 h of reaction time, which implied less influence of catalyst loading on CMEC yield between the reaction time of 2 to 6 h. However, a maximum CMEC yield of 68% was achieved at higher catalyst loading and reaction time of 12% (w/w) and 12 h, respectively. A declining effect was observed in Figure 14a as the catalyst loading goes above 12% (w/w). This reflects that the optimum catalyst loading had been exceeded. A similar trend was reported by Onyenkeadi and research team on declining effect of catalyst loading beyond the optimum reaction time. Increase in the amount of catalyst loading can increase the number of active sites on the catalyst surface, and consequently, increases number of radicals (see supplementary information sheet, Figure S1). However, excessive increase of catalyst concentration beyond the optimum reaction time can result in a catalyst deactivation. This phenomenon is totally in agreement with the recent reports of Feilizadeh et al. [69].
Energies 2020, 13, x FOR PEER REVIEW 19 of 27 contour plots show less curvature up to 7 h of reaction time, which implied less influence of catalyst loading on CMEC yield between the reaction time of 2 to 6 h. However, a maximum CMEC yield of 68% was achieved at higher catalyst loading and reaction time of 12% (w/w) and 12 h, respectively. A declining effect was observed in Figure 14a as the catalyst loading goes above 12% (w/w). This reflects that the optimum catalyst loading had been exceeded. A similar trend was reported by Onyenkeadi and research team on declining effect of catalyst loading beyond the optimum reaction time. Increase in the amount of catalyst loading can increase the number of active sites on the catalyst surface, and consequently, increases number of radicals (see supplementary information sheet, Figure S1). However, excessive increase of catalyst concentration beyond the optimum reaction time can result in a catalyst deactivation. This phenomenon is totally in agreement with the recent reports of Feilizadeh et al. [69].

Interactive Effect of Catalyst Loading and Pressure
The exponential interaction effect between catalyst loading and pressure at a constant reaction time of 12 h and a temperature of 353 K is presented in Figure 15. However, the interaction produced a different effect on CMEC yield at different levels of interaction (i.e., different levels of interaction produce different effects on the ECH conversion). For example, Figure 15 shows that at the start of the reaction, 5% (w/w) of catalyst loading at 7 bar of CO 2 pressure produced an increasing effect on the CMEC yield. As the catalyst loading was further increased from 5% to 10% (w/w), the CMEC yield was observed to increase steadily from 40% to 68% corresponding to an increase in CO 2 pressure from 7 to 11 bar. The CMEC yield was highest (68%) at a maximum catalyst loading of 12 % (w/w), when the CO 2 pressure was maintained at 11 bar as shown in Figure 15. However, a negative effect of excessive increase in CO 2 pressure was observed on CMEC yield (a drop to 64%) at this level of interaction between catalyst loading and pressure. This phenomenon can be attributed to catalyst deactivation at increased CO 2 pressure beyond the optimum. A similar experience was reported earlier by Zhang et al. [70]. The group have recorded a higher propylene carbonate (PC) yield with a fixed amount of immobilized ionic liquid/ZnCl 2 at a CO 2 pressure of 1.5 MPa, however, a lower PC yield was observed at a higher CO 2 pressure of 2 MPa. Furthermore, they claimed that this phenomenon occurs when acidic CO 2 dissolves in basic epoxide to form a liquefied CO 2 −epoxide complex, thereby inducing catalyst deactivation. Energies 2020, 13

Multiobjective Process Optimization
The growing quest for greener substitute for fossil fuel has led to increased production, process optimization, and application of organic carbonate. As a result, the use of RSM has received more attention over conventional optimization methods in order to investigate process optimum

Multiobjective Process Optimization
The growing quest for greener substitute for fossil fuel has led to increased production, process optimization, and application of organic carbonate. As a result, the use of RSM has received more attention over conventional optimization methods in order to investigate process optimum conditions and the interactive relationships between effective working variables. Although, finding the optimal reaction parameters for a single response using RSM is relatively simple; however, the optimization of several responses at the same time is not an easy matter. Table 5 shows the optimization targets for this study, the values have been set to maximize the process productivity. Targets for both ECH conversion and CMEC yield have been set to reach the maximum values while both the reaction temperature and time have been targeted to minimum values with a viewpoint of reducing production cost at a maximum economic gain. Because of the catalyst efficiency and stability at optimum conditions, as a result, no specific target has been set for catalyst loading. Based on the models generated and the accuracy between the actual experimental and predicted results, it can be construed that the model shows high consistencies between the two results where the relative errors of the predicted results from the experimental data are 1.55% and 1.54% for ECH conversion and CMEC yield, respectively. The similarity between the predicted and experimental results at the optimum conditions has validated the predicted optimum conditions. The experimental results concluded that increase in reaction parameters increases ECH conversion and CMEC yield being 93% and 68%, respectively.

Catalyst Reusability Studies
In view of large scale industrial applications and to minimize production cost, the reusability studies of Zr/ZIF-8 catalyst have been investigated. The catalyst reusability process has also followed strict eco-regulation after all the predicted optimum parameters have been derived from BBD of RSM. The experiments were carried out in a high-pressure reactor at optimum reaction conditions, at 353 K, 11 bar with fresh 12% (w/w) ZIF-8 catalyst loading, for 12 h, and at a stirring speed of 350 rpm. The catalyst after Run 1 in the cycloaddition reaction was washed with ethanol and acetone, centrifuged, and oven dried at 343 K for 12 h before reuse. The recovered catalysts were reused for up to seven subsequent experiments as shown in Figure 16, following the same experimental procedure. The catalyst exhibited no loss of activity indicating the catalyst stability for cycloaddition reaction of CO 2 epichlorohydrin. Incorporating zirconium into ZIF-8 has significantly increased the catalytic performance of Zr/ZIF-8 with the conversion of ECH and the yield of CMEC being 93% and 68%, respectively. The activity of reused Zr/ZIF-8 catalyst showed consistent stability over seven subsequent runs as indicated in Figure 16. Although, a very slight decrease in the yield of CMEC from 68% (fresh) to 67% (recycled) was observed in the seventh run. Carbonaceous material formed during the reaction may explain in part the lower activity of the recycled catalysts [71].

Conclusions
In this study, Zr/ZIF-8 catalyst has been successfully used for process optimization in the synthesis of CMEC using RSM. In total, 29 runs of experiments were conducted for optimum design and modelling. The developed model was validated to assess the agreement between its predictions and a set of experimental data. The development of a novel Zr/ZIF-8 catalyst via a simple low cost solvothermal method has demonstrated that the catalyst is viable for large-scale industrial applications. The catalyst has shown a good substrate tolerance as demonstrated by its activity towards epichlorohydrin. More importantly, the reaction has been carried out under solvent free and co-catalyst free conditions. The heterogeneity of the catalyst has been proven by recovering and reusing the catalyst for up to seven times without any significant loss in catalytic activity. Furthermore, PXRD, FT-IR, and TGA analysis (see supplementary information sheets, S2, S3 and S4) of the recycled catalyst shows that the catalyst framework is quite stable after recycled experiments. The high selectivity towards epichlorohydrin carbonate, simple separation of catalyst by centrifugation, and excellent recyclability demonstrated that the catalyst is viable for industrial applications. We believe that this work could provide a new direction for designing more sustainable heterogeneous catalysts for greener synthesis of organic carbonates via CO2 utilization.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: Proposed free radical mechanism for cycloaddition reaction of CO2 and ECH, Figure S2: X-ray diffraction (XRD) patterns of (a) ZIF-8, (b) Zr/ZIF-8 and (c) recycled Zr/ZIF-8 catalysts, Figure S3: FTIR spectra of (a) ZIF-8, (b) Zr/ZIF-8 and (c) recycled Zr/ZIF-8 particles, Figure S4 Although, the difference in the error bars status between the ECH conversion and CMEC yield may be statistically significant, this may be attributed to the formation of some side products associated with the coupling reaction of CO 2 and ECH. The following side products have been identified by the GC analysis; 3-chloropropane 1,2-diol and 2,5-bis (chloromethyl)-1,4-dioxane.

Conclusions
In this study, Zr/ZIF-8 catalyst has been successfully used for process optimization in the synthesis of CMEC using RSM. In total, 29 runs of experiments were conducted for optimum design and modelling. The developed model was validated to assess the agreement between its predictions and a set of experimental data. The development of a novel Zr/ZIF-8 catalyst via a simple low cost solvothermal method has demonstrated that the catalyst is viable for large-scale industrial applications. The catalyst has shown a good substrate tolerance as demonstrated by its activity towards epichlorohydrin. More importantly, the reaction has been carried out under solvent free and co-catalyst free conditions. The heterogeneity of the catalyst has been proven by recovering and reusing the catalyst for up to seven times without any significant loss in catalytic activity. Furthermore, PXRD, FT-IR, and TGA analysis (see supplementary information sheets, S2, S3 and S4) of the recycled catalyst shows that the catalyst framework is quite stable after recycled experiments. The high selectivity towards epichlorohydrin carbonate, simple separation of catalyst by centrifugation, and excellent recyclability demonstrated that the catalyst is viable for industrial applications. We believe that this work could provide a new direction for designing more sustainable heterogeneous catalysts for greener synthesis of organic carbonates via CO 2 utilization.