Croton gratissimus oil was produced in the laboratories at Mangosuthu University of Technology (MUT) in an extraction process of Croton gratissimus grains (Figure 1) carried out using n-hexane as an extraction solvent (Jiyane et al., 2018). The oil produced had an acid value of 21.46 mg KOH/g oil and a specific gravity of 0.928 (Table 1). Sulfated zirconia catalyst, SO₄ ^2–/ZrO₂, was synthesized in the laboratories at MUT following the method proposed by Chen et al. (1993) and modified by Reddy et al. (2005). Potassium hydroxide (KOH ≥ 85%) and zirconyl chloride octahydrate (ZrOCl₂·8H₂O), a precursor for the sulfated zirconia catalyst SO₄ ^2–/ZrO₂, were purchased from Sigma-Aldrich, South Africa. Methanol AR of 99% purity was purchased from Merck Laboratories in South Africa.

An analytical balance, RADWAG® AS 220/C/2, was used for accurate weight measurement of catalyst samples used in FAME production. A HERMLE centrifuge was used for final separation of the esterification and transesterification products. The laboratory-assembled equipment, shown in Figure 2, was used in the esterification and transesterification reactions. A PerkinElmer Clarus® 580 Gas Chromatograph equipped with a Clarus® SQ 8S Mass Spectrometer (GC-MS) was used in the characterization of finished biodiesel product.

Biodiesel is mainly produced in a transesterification reaction, a chemical reaction where alcohol from an ester is displaced by another alcohol in a process similar to hydrolysis, except that transesterification uses an alcohol instead of water. During this reaction, triglycerides from oils and fats react with an alcohol in the presence of a catalyst to produce fatty acid esters (biodiesel) and the by-product glycerol. It is important to note that with recent advances in research in this field, other methods have been developed that do not require the use of a catalyst, that is, catalyst-free. In this regard, supercritical methods that employ supercritical methanol have been successfully used in biodiesel processing (Madras et al., 2004) albeit being energy intensive and inherently unsafe. In a catalytic environment, a basic or acidic catalyst is used to convert the glycerol-based esters (triglycerides), which make up fats or oils, to alcohol-based esters (methyl esters) giving off free glycerol as a by-product. Glycerol, after separation from the reaction mixture, may be purified and sold for use to the pharmaceutical, cosmetics, foods, and the plastics industries (Vicente et al., 1998). The mass balance of croton gratissimus production can be written as follows:

1,000 kg + 110.50 kg → 1,004.64 kg + 105.86 kg

Fat oil + 3 Methanol → 3 Methyl Ester + Glycerol

This mass balance is based on the fatty acid composition of Croton gratissimus oil shown in Table 2 and assumes complete conversion of triglycerides to fatty acid methyl esters. It therefore gives the theoretically attainable FAME yield from this virgin oil.

The calcined catalyst was characterized using the X-ray diffraction (XRD) analyses to establish the phases and measure the crystalline sizes; the scanning electron microscopy (SEM) to evaluate the catalyst’s surface topography and composition, and the transmission electron microscopy (TEM) that provided morphological, compositional, and crystallographic information on the catalyst.

Supplementary Figure 1 provided in supplemental materials shows the X-ray diffraction profile obtained from the measurements carried on the SO₄ ^2–/ZrO₂ catalyst. The main phase was found to be ZrO₂, which matched the tetragonal lattice pattern (red lines). The tetragonal phase was predominant at 2θ of 30.5°, and Baddeleyite, a monoclinic ZrO₂ phase (powder blue lines), was mostly found at 2θ of 28.2°. The presence of both phases (with a dominant monoclinic phase) in the crystal lattice indicated that the doped sulfate ions did not successfully influence the phase modification of zirconia from the thermodynamically more stable monoclinic phase to the metastable tetragonal phase. For this reason, the synthesized SO₄ ^2–/ZrO₂ catalyst was observed to be highly active in catalyzing only the esterification reaction and not the transesterification reaction.

Morphological studies on the catalyst were performed using the TEM and its composition was analyzed using the SEM. Supplementary Figure 2 (Supplemental materials) shows the morphology of SO₄ ^2–/ZrO₂ catalyst calcined at 620°C for 4 h. Highly packed crystal morphology with irregularly arranged ZrO₂ aggregation was observed. There was a high loading amount of Zr in the catalyst structure as proven by a 61.78 mass% Zr from the elemental analysis of the catalyst (Table 3). The presence of all the constituent atoms (Zr, O, and S) that make sulfated zirconium catalyst was confirmed by spectrum 1 in Supplementary Figure 3 (Supplemental materials). The spectrum was developed from a randomly picked area on the face of the sample (Supplementary Figure 3). Because of its rather low mass and abundance in the catalyst, the concentration of sulfur was low and in some spot spectra (spectra 2 and 5 on Table 3) of Supplementary Figure 3 remain undetected (Figure 6). To determine the concentration of sulfates and, hence, the acidity and activity of sulfated zirconium catalyst, other methods are usually employed. The mean diameter of the monoclinic zirconia crystals is ± 30 nm while that of the tetragonal crystals is ± 200 nm (Miranda M et al., 2015). In Supplementary Figure 4, the TEM micrographs show the sizes of most of the crystals in the SO₄ ^2–/ZrO₂ sample to range between 10 and 100 nm, indicating an overwhelming presence of the monoclinic phase.

A PerkinElmer Clarus® 580 Gas Chromatograph equipped with a Clarus® SQ 8S Mass Spectrometer (GC-MS) was used to measure the percentage purity (% P) of FAME samples. The gas chromatograph (GC) used has a 30 m length PerkinElmer Elite-5MS II capillary column of 0.5 µm film thickness and 0.25 mm inner diameter. The stationary phase in the column selected is a low polarity 5% diphenyl/95% dimethyl polysiloxane, and the mobile phase is pure helium (He), flowing at 1 ml/min, used as carrier gas.

The FAME sample of ±100 µL was dissolved in 10 ml of ethanol and a 1 µL volume of the diluted sample was injected into the GC. Methyl heptadecanoate was used as an internal FAME standard. The oven temperature was initially kept at 140°C for 5 min and ramped up at 4 °C/min to a maximum of 240°C for a total run time of 31 min. The sampling rate was maintained at 1.5625 pts/s enabling the detector to capture enough data points across very narrow peaks. The purity of FAME was then determined from the chromatographic peak areas using the following formula (Jeong et al., 2009): Percentage  FAME Purity ( % P ) = ( ∑ A FS ) − A IS A IS × C IS × V IS m FS × 100 , (3) where, ∑ A FS is the sum of all the methyl ester (FAME) peak areas, A IS is the peak area of the internal FAME standard (methyl heptadecanoate), C IS is the concentration in mg/ml of the solution of the internal standard, V IS is the volume in ml of the internal standard solution, and m FS is the mass in mg of the FAME sample.

The FAME analysis method described above, that employs Eq. 3, was adopted in this study as it is the most preferred when the quantity (purity) rather than the quality (composition) of FAME is being investigated in biodiesel synthesis. But, when the composition of FAME in the biodiesel test sample is of an essence, FAME standards of known concentration are used to plot calibration curves which are then used in the determination of individual methyl esters.

The concentration and volume of the internal standard (methyl heptadecanoate) were 33 mg/ml and 100 μL, respectively. The peak area ( A IS ) of methyl heptadecanoate obtained was 98,147 × 10³ units. Supplementary Figure 5 (Supplemental materials) represents the GC-MS spectrum used for the analysis of FAME purity in the 107.50 mg biodiesel sample injected. The fatty acid methyl esters in the sample were evaluated using the GC-MS library that identified the most prominent methyl ester on each peak. The presence of the following methyl esters viz., oleate, palmitate, stearate, linoleate, palmitoleate, arachidate, and linoleneate and their peak areas (Supplementary Figure 5) were noted on each spectrum. Percentage FAME purity was then calculated from the total peak areas (Table 4) using Eq. 3 (Jeong et al., 2009).

At the end of the experiments, the data collected was subjected to an optimization study to find the optimum operating conditions for both the esterification and the transesterification processes. Henceforth, the optimization studies of these processes will be referred to as OPTIMA1 and OPTIMA2, respectively. An optimization study further carried out to find a set of optimum independent variables that could simultaneously satisfy both optimum values of FAME yield and FAME purity obtained in OPTIMA2 will be referred to as OPTIMA3. Thereafter, the response surface methodology (RSM), that employs the central composite design (CCD) and analysis of variance (ANOVA) techniques, was employed in these studies.

The CCD was applied with three design factors, viz., catalyst concentration ( X 1 ) , the methanol-to-oil ratio ( X 2 ) , and the reaction temperature ( X 3 ) . These factors were selected to optimize the chosen responses of acid value ( Y A ) in OPTIMA1, and FAME yield ( Y Y ) and FAME purity ( Y P ) in OPTIMA2. The selection of levels (uncoded variables) for each factor was based on the literature related to the optimization of similar biodiesel processes. Processes of interest were those that draw their feedstock from high FFA non-edible vegetable crops and use the two-step production process approach, that is, esterification and transesterification reactions.

For both processes (OPTIMA1 and OPTIMA2), 2³ full-factorial rotatable CCD designs for three independent variables at five levels were employed and the total number of experiments conducted was 2 × 20, including axial and replicate experimental runs. Experiments were performed in a randomized order to balance out the effect of any nuisance variable that may influence the observed response (Montgomery and Runger, 2010). The Design-Expert® version 10 software was used to design the experiments for regression and graphical analyses of the data obtained. Acid values obtained in OPTIMA1 and the FAME yield and purity obtained in OPTIMA2 and OPTIMA3 were all taken as responses of the designed experiments.

Statistical analysis of the model was performed to evaluate the analysis of variance (ANOVA). The Fisher’s F-test was performed to test the significance of adding quadratic terms to the two-factor interaction (2FI) model. Ideally a small p-value (Prob>F less than 0.05) is desirable as it indicates that the addition of these (quadratic) terms has contributed to an improvement in the model. The goodness of fit of the model was evaluated by the determination coefficient (R ²), adjusted determination coefficient (Adj. R ²), and the coefficient of variation (C.V.).

The coefficient of variation (C.V.) or most commonly known as the relative standard deviation is the standard deviation expressed as a percentage of the mean and is crucial in measuring the reproducibility of the model. Predicted R ² is a measure of how well the model predicts a response value and adjusted R ², on the other hand, measures the amount of variation around the mean that is explained by the model. For the model to fit the collected data, the predicted R ² (Pred. R ²) value must be in reasonable agreement with the adjusted R ² (Adj. R ²) value by not more than 0.2 point difference. The adequate precision (Adeq. Precision) measures the signal-to-noise ratio and should ideally be greater than 4. A value greater than 4 is desirable as it indicates that the developed model has a high degree of precision and can thus be used to navigate the entire design space.

From the developed quadratic regression model, a combination of independent variables ( X 1 , X 2 , and X 3 ) , forming an optimal set, that fell within the region of 95% confidence interval (95% CI) of the chosen response ( Y A , Y Y , or Y P ) was taken as the optimum set of variables that is most likely to produce that optimum response. To narrow down the search for the best compromise within the optimum region, comprising a multitude of independent variables matching multiple responses, optimum acid value in OPTIMA1 and FAME yield and FAME purity in OPTIMA2 were targeted and one independent variable set as a constraint.

In OPTIMA1, an acid value of 2.963 and a reaction temperature of 64°C were used as constraints to reduce the number of plausible optimum set. FAME yield of 84.51% and FAME purity of 90.66% obtained in OPTIMA2 and reaction temperature of 63.5°C were selected as constraints for use in OPTIMA3. The choice of the final optimum set of independent variables was based on process economics, that is, reduced temperatures and low levels of catalyst usage.

Overlay plots were drawn, in graphical optimization, to display the areas of feasible response values within the factor space, and with “flag inserts” to show the optimum solution set. In the overlay plots shown in Figure 3, the region colored in green (1) represents a region where graphical optimization solutions that meet the specific criteria can be found. The gray region (2) represents the graph area with solutions outside of model prediction. The region colored in red (3) represents the graph area outside the interval bound for at least one response.

High acid value of Croton gratissimus oil (21.455 mg KOH/g) necessitated that an esterification reaction be conducted prior to the main transesterification reaction. This was done in order to reduce the acid value of the extracted oil to below 4 mg KOH/g of oil. The effects of catalyst concentration ( X 1 ) , the methanol-to-oil ratio ( X 2 ) , and the reaction temperature ( X 3 ) on the acid value were investigated. They were each varied in 5 levels as shown in Table 5. The values of the levels were established through the OFAT esterification experiments conducted over the synthesized monoclinic SO₄ ^2–/ZrO₂ catalyst.

Table 6 gives the results of acid value reduction for the 20 experimental runs carried out to find the optimum conditions for the esterification process. To evaluate the accuracy of the model in predicting the response values that match the actual experimental values, various statistical parameters are provided in Table 7.

The developed model was able to predict most of the response values as the points on the plot lie in close proximity to the 45° diagonal line (Figure 4).

The experimental data obtained in the esterification process was analyzed by the ANOVA and is best described by the following quadratic regression model: Y A = 17.362 − 0.456 X 1 − 0.152 X 2 − 0.239 X 3 − 2.958 × 10 − 4 X 1 X 2 − 2.854 × 10 − 3 X 1 X 3 + 9.917 × 10 − 4 X 2 X 3 + 0.020 X 1 2 + 1.254 × 10 − 3 X 2 2 + 1.670 × 10 − 3 X 3 2 , (4) where Y A is the acid value (mg KOH/g) of the esterified oil and X 1 , X 2 , and X 3 are uncoded values of the independent variables, viz., catalyst concentration (mass%), the methanol-to-oil ratio, and the reaction temperature (°C), respectively. An increase in either catalyst concentration, the methanol-to-oil ratio, or the reaction temperature resulted in a reduction in the acid value of croton gratissimus oil. The negatively correlated linear terms in the developed quadratic regression model (Eq. 4) are all evidence to this assertion. The quadratic model also revealed that catalyst concentration was the most significant independent variable (p-value < 0.0001), exhibiting a larger negative value to its linear coefficient term than other model factors terms.

The ANOVA results in Table 7 showed the significant nature of the developed quadratic regression model with a reasonably low probability value, Prob > F: (p-value) < 0.0001. The Fisher’s F-test gave an F value for the model of 55.94, indicating a less than 0.01% chance that an F value this large could be due to noise. The lack-of-fit with an F-value of 1.23 implied that the lack-of-fit was not significant and was proof that the model developed adequately represented the collected data. The value of the coefficient of determination (R ²) was found to be 0.9805, indicating that only 0.20% of the total variations was not explained by the developed regression model. The predicted R ² of 0.9045 was in reasonable agreement with the adjusted R ² of 0.9630 with only an acceptable difference of less than 0.20. The developed model had a high degree of precision as the adequate precision found was high enough at 28.141. The coefficient of variation (C.V.) of 3.70% indicated that the model was reproducible.

From the developed model, an acid value of 2.693 mg KOH/g could be obtained when operating the process under optimum conditions of 10.96 mass% catalyst concentration, 27.60 methanol-to-oil ratio, and a reaction temperature of 64°C. The 91% reduction in the acid value of Croton gratissimus oil (based on the lowest experimental value obtained) achieved in this esterification step could positively influence the reaction extent in the next transesterification process step. Figure 4 also shows surface plots of the effects of catalyst concentration, the methanol-to-oil ratio, and the reaction temperature on the acid value of Croton gratissimus oil during the esterification process. The amount of catalyst used had the most significant effect on the acid value of oil. This is shown by the almost flat-shaped response surfaces obtained for the methanol- to-oil ratio and reaction temperature on the acid value when the catalyst concentration was kept constant at 10.96 mass% (Figure 4). The focus of the esterification reaction was on the use of the sulfated zirconia, SO₄ ^2–/ZrO₂ in its current form as a predominantly monoclinic phase and evaluating its effectiveness in converting free fatty acids to di and monoglycerides. For this reason, a much wider range of catalyst concentration, between 6 and 14 mass%, was selected for the study whilst feeding methanol in excess (more than the stoichiometric amount). Arora et al. (2015) in their study on the esterification of high free fatty acid rice bran oil found out that there was no significant increase in the FFA conversion when the methanol-to-oil ratio was increased beyond 20. In the current study, methanol was fed in a range of between 20 and 40:1 to rule it out as the contributing factor, should there be low or no conversion of FFA at the end of the reaction.

The amount of the SO₄ ^2–/ZrO₂ catalyst used in the esterification reaction was found to have the most significant effect in the reduction of the acid value of Croton gratissimus oil from 21.455 mgKOH/g to an average of 2.927 mg KOH/g oil. The highest and lowest acid values were obtained when the same operating conditions of temperature (70°C) and the methanol-to-oil ratio (32.50) whilst more than doubling the catalyst concentration from 6 to 14 mass% (Table 6 standard experimental run 9 and 10). This dependency of the esterification reaction on the catalyst concentration was consistent with the observations by Kumar Tiwari et al. (2007). Their optimized RSM model for biodiesel produced from high FFA Jatropha curcas oil showed that catalyst concentration was the most significant factor (followed by the methanol-to-oil ratio) in the conversion efficiency of free fatty acids (FFA) in oil.

Within the lower limit of the 95% confidence interval, the developed quadratic model predicted an optimum acid value of 2.693 mg KOH/g of oil when the esterification process was operated at 10.96 mass% catalyst concentration, 27.60 methanol-to-oil ratio, and 64°C reaction temperature. In a two-step biodiesel production process, Muthu et al. (2010) converted Azadirachta indica A. Juss oil to biodiesel using SO₄ ^2–/ZrO₂ and KOH as catalysts in the esterification and transesterification reactions, respectively. From an unoptimized process operating at the methanol-to-oil ratio of 9:1, catalyst concentration of 1%, reaction time of 2 h, and reaction temperature of 65°C, they obtained a 94% FFA conversion efficiency in the esterification step with a predominantly tetragonal-phased SO₄ ^2–/ZrO₂ catalyst.

Independent variables selected for an investigation in the optimization of the transesterification process were catalyst concentration ( X 1 ) , the methanol-to-oil ratio ( X 2 ) , and the reaction temperature ( X 3 ) . They were each varied in 5 levels, that is, − α ,   − 1 ,   0 ,   + 1 ,   + α , as shown in Table 5, to establish their effects on the responses of FAME yield and FAME purity. A response surface methodology’s 2 ^k full-factorial central composite rotatable design (CCRD) was employed giving a total of 20 experiments, evaluated from 2 ^k + 2k + n _c , where k = 3 and n c = 6. Experiments were carried out in a randomized order to minimize the effects of uncontrolled factors. The final experimental results are presented in Table 8. To evaluate the analysis of variances (ANOVA), statistical analyses of the 2 developed quadratic models were carried out. Optimum operating conditions that gave the highest percentage FAME yield (Response 1) and FAME purity (Response 2) were evaluated.

The two quadratic models developed for the transesterification process predicted two sets of factor values for catalyst concentration, the methanol-to-oil ratio, and reaction temperature corresponding to the optimum values of the 2 responses, FAME yield and FAME purity. Because the transesterification process is but one process, it must be described by one set of variable values for the attainment of these optimum responses. A combined model was, therefore, developed for the optimum set of independent variables that gave responses of FAME yield and FAME purity within the region of 95% confidence interval (95% CI). This set of independent variables had to simultaneously satisfy the optimum responses obtained in OPTIMA2.

Numerical optimization techniques, that use the response targeting approach bound by constraints and the graphical optimization technique that utilizes overlay plots to display regions of feasible responses that fit the optimization criteria, were employed. This task was achieved by overlaying critical response contours on a contour plot (overlay plot) after targeting optimum responses from OPTIMA2 and using constraints to narrow down the search for the optimum set. Reaction temperatures ranging between 63.50 and 64°C (below 64.50°C, the boiling point temperature of methanol) and low catalyst concentration between 10.50 and 10.98 mass% for OPTIMA1 and 1.40 and 1.60 mass% for OPTIMA3, were used as constraints. The choice of these constraints was based on process economics. From the combined model optimization, it was found that when the transesterification process is operated at 1.439 mass% catalyst concentration, 7.472 methanol-to-oil ratio, and 63.50°C reaction temperature, 84.51% FAME yield and 90.66% FAME purity could be achieved. These independent variables were, therefore, taken as the optimum set that satisfies optimum responses predicted by the developed models.