Influence of moisture content, temperature, and time on free fatty acid in stored crude palm oil

Consequent to the importance of crude palm oil (CPO) to global food processing industries, and the need for quality assurance of CPO. A kinetic model that describes changes of free fatty acid (FFA) in industrially stored CPO has been developed. CPO FFA is a well-known indicator of the deterioration of CPO. The effect of initial moisture content, storage temperature, and time on CPO FFA have been investigated in this work. Specifically, statistical multi-regression models for changes in FFA and moisture content (MC) were developed at P-value < 0.05 or 95% confidence interval fence. It was found that CPO FFA increases with an increase in moisture content, temperature, and time in their linear term and in respect to decreases in their quadratic term, and interaction between moisture content and temperature. The CPO MC was also found to decrease with an increase in temperature and time and increases in the quadratic term of temperature. Although while the model for CPO FFA, based on Fisher's F-test: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{F}}_{\mathrm{model}}(6.80)<{\mathrm{F}}_{95\mathrm{\%}}(19.30)$$\end{document}Fmodel(6.80)<F95%(19.30), showed no lack-of-fit; that of CPO MC showed lack-of-fit, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{F}}_{\mathrm{model}}(13.67)\nless {\mathrm{F}}_{95\mathrm{\%}}(4.39)$$\end{document}Fmodel(13.67)≮F95%(4.39). Furthermore, based on inference from the statistical model, their kinetic models were also developed. While the CPO FFA kinetic, found to be a half-order kinetic model and its other auxiliary models showed a very good fit (R2 {0.9933–0.8614} and RMSE {0.0020–3.6716}); that of CPO MC was a poorly fitted first-order kinetic model (R2 {0.9885–0.3935} and RMSE {0.0605–17.8501}).

www.nature.com/scientificreports/ that there are more triglycerides with palmitic chains than other triglyceride constituents, and as such palmitic acid significantly influences CPO FFA. Almeida et al. 17 , and Taluri et al. 18 also reported that storage temperature influenced the FFA of CPO and olive oil respectively. Furthermore, Lin et al. 19 reported a kinetic model for FFA formation in lipid extracted from stored almonds, with consideration of the influence of relative humidity, and temperature. Although having highlighted from literature that the MC and temperature of stored vegetable oil increase FFA value, there are however no reports on a kinetic model to predict CPO FFA changes based on these two highlighted storage factors. This model when developed would facilitate the simulation, monitoring, and control of FFA in industrially stored CPO. The dynamics of FFA in vegetable oil can be investigated through the reactants (i.e. reaction of glyceride and water molecules), Eq. (1) [20][21][22][23][24] . However, the measurements in this approach require the simultaneous analysis of glycerides, water, and FFAs, using expensive and complex equipment such as Gas Chromatography, Liquid Chromatography, and HPLC [25][26][27] . Therefore a simplified approach based on the product (i.e. change in CPO FFA, % FFA via titration method), Eq. (2) is considered in this work. Furthermore, the model development approach for Eq. (1) is suggestive to be more complex than Eq. (2), as it can constituent multiple subequations of the hydrolysis process 20 , and as such more variables. Where, r FFA is the reaction kinetic, y FFA = % FFA , k is the reaction constant, n is the reaction order, i is the specific glyceride molecule considered, x gly,i and x water are the compositions of glycerides, and water in the CPO.
Therefore, in this work, CPO FFA and MC with known initial FFA and MC stored at different temperatures in a surrounding of known relative humidity will be analysed at different storage times. The objectives shall involve: Analysis and collection of adequate quantity of CPO of specific MC from commercial CPO processing company; Design of experimental via Box-Behnken design (BBD) of the minimum and maximum limit of factors (i.e. initial MC (with known measured FFA), temperature and time) and with final FFA and MC as responses; Storage of the CPO samples at the designed BBD condition, and subsequent evaluation of the specified responses; Deduce statistical multi-regression model for the designed BBD, and evaluation of the factors and equations at 95% confidence level using the Fisher's F-test; Evaluation of FFA and MC kinetics and deduction of an n th order kinetic model.

Materials and methods
Collection of the source sample. CPO samples were selected from the regular sampling routine of a CPO processing company based on the acceptable maximum limit (0.20%) and range (0.015-0.30%) of MC in edible vegetable oil 14 . Three sample sets of 0.20 ± 0.02%, 0.25 ± 0.02%, and 0.30 ± 0.02% MC were identified, and their FFA were also determined (5.17, 5.05, and 4.28% respectively). Thereafter adequate quantity (1000 ml) of these samples were collected, sealed, and stored as a source sample in an enclosure with an Ultraviolet-C radiation lamp to limit FFA increase as reported by Said et al. 28 . Table 1 shows the properties of the CPO source samples given by the company.
Determination of free fatty acid. Based

Statistical model analysis.
Design of experiment via BBD was implemented to determine whether MC, temperature, and storage time indeed influence the CPO FFA as reported in literature. This was achieved by the deduction of multi-variable regression model using Python Jupyter Notebook, and was also deduced for CPO MC.

Design of experiment.
A "three-variable-three-level" BBD of three (3) centre points was adopted resulting in fifteen (15) experimental runs, N , with the three specified variables considered for each response. The design analyses the contribution of these variables in terms of linear, quadratic, and interaction effects in the prediction of CPO FFA and MC via a generalised second-order polynomial regression model, Eq. (9). Equation (8) and Table 2 shows the relationship between coded, x and actual values, ξ of the specified variables for the analysis. Where ξ = 0.5 (ξ max + ξ min ) , ξ max and ξ min are the mean, maximum and minimum actual value of the experiment.
In a similar manner for which the three sample sets for CPO have been chosen (i.e., based on the CPO MC standard for edible vegetable oil). The upper and lower limit of temperature has also been chosen based on the suggested range at which CPO is processed and stored by commercial CPO processing companies.
Where Y , is the response, i.e. % FFA or % MC as given by Eqs. (5) and (7) respectively. x i and x j represent the specified variables of consideration, β 0 is the constant of the model, β i is the linear term coefficient,β ii is the quadratic term coefficient, β ij is the interaction coefficient and n is the number of variables considered.
Reaction kinetics analysis. The reaction kinetics for the CPO FFA and MC was deduced using Eq. (10), by curve fitting experimental data of FFA and MC for the three source samples collected. 50 g of each of these source samples was poured into a 50 ml beaker of height, 60.96 mm and diameter, 38.10 mm. The sample was left opened and placed in an oven, and the analysis of the sample's FFA and MC following the procedures described earlier were performed at a temperature interval of 10 ℃ from 35 to 85 ℃, and at a time interval of 6 h from 6 to 120 h. The rate constant, k for the reaction kinetics were deduced on the assumption that the reaction order (i.e.,n , the power factor on % FFA , y FFA and % MC , y MC ) is within zero-to-second order (i.e., n = 0, 0.5, 1.0, and 2.0).
In modelling FFA kinetics, r FFA , it is assumed k is a function of temperature, and change in initial moisture content, MC i , such that it is a product of the Arrhenius relationship based on temperature and MC i (Where MC i = MC i − MC i,0.20% i.e. changes in MC of a sourced sample from the standard 0.20% MC) as expressed by Eqs. (11) and (12).  (11) is deduced by curve fitting the value of k computed for each CPO source sample at different temperatures based on the results from Eq. (10). While k m , Eq. (12) is deduced by curve fitting the value of the ratio of k 0 in Eq. (11) with respect to k 0 of the maximum allowable 0.20% MC, i.e. k m = k 0 /k 0,0.20% . This manipulation implies reexpressing k as Eq. (13). Where R is the ideal gas constant, 8.314 J·mol −1 ·K −1 , T(K) is the temperature, E(J·mol −1 ) is the activation energy, ϕ(J·mol −1 ) is activation energy due to MC i , and 273 is the temperature in Kelvin at standard temperature and pressure.
In modelling MC kinetic, r MC , Eqs. (10) and (11) are the only applicable equations to utilise.
Evaluation of deduced kinetic model. The goodness-of-fit for the kinetic models deduced were evaluated using R-squared (R 2 ), Eq. (14) and Root mean squared error (RMSE), Eq. (15). The use of these two evaluation criteria is hinged on the report that R 2 is inappropriate when used for demonstrating the performance or validity of certain nonlinear models 33 . Consequent to this report Jim 34 has suggested that RMSE would be appropriate for such cases. On the basis that the kinetic models would be evaluated for n = 0, 0.5, 1, and 2, hence the possibility of nonlinearity. Hence, for nonlinear models, RMSE would be used as the main criterion even though R 2 would also be reported. The R 2 value ranges between 0 and 1, the closer to unity the better the model, while RMSE ranges between 0 and ∞ , the closer to zero the better. Where y and y is the output and mean output of the experimental data, y is the curve fitted model output, N is the number of the experimental run, and p is the number of dependent variables in the model.
The curve fitting of kinetic models was implemented using Lsqcurvefit curve-fitting tools in MATLAB.

Results and discussion
Result of design of experiment. Before developing stored CPO FFA and MC kinetic models, statistical analysis via the Box-Behnken design of experiment was used to investigate the variables that significantly influence them. This analysis was achieved by comparison of coefficients, β 0 , β i , β ii and β ij computed with P-value ** < 0.05 or 95% confidence interval fence. The resulting multi-regression models for the statistical analysis are given by Eqs. (16) and (17). The model shows, FFA increases with an increase in initial moisture content, x 1 , temperature, x 2 , and time, x 3 in their linear term and decreases in their quadratic term, and interaction between moisture content and temperature, x 1 x 2 . This result is quite similar to Gawrysiak-Witulska et al. 35 reported on the effect of moisture content and temperature on the degradation of phytosterol, which is also found in lipid. In the Gawrysiak-Witulska et al. 35 report a three factorial analysis of variance of the factors given in Table 2 was also performed, but only in respect to linear ( x 1 , x 2 , and x 3 ), and interactive terms ( x 1 x 2 , x 1 x 3 , x 2 x 3 and x 1 x 2 x 3 ), and it was found that moisture content, temperature, storage time significantly influenced the degradation of rapeseed phytosterol in all highlighted terms.
Furthermore, CPO MC in reference to Eq. (17) decreased with an increase in the linear terms of temperature, x 2 , and time,x 3 , and increased with an increase in the quadratic term of temperature, x 2 2 .
Although these equations showed the significant variables and their relationship to FFA and MC. While y FFA showed no lack-of-fit, i.e. f model (6.80) < f 95% (19.30) , y MC showed lack-of-fit, i.e. f model (13.67) ≮ f 95% (4.39) , as such y FFA is adequate for prediction of FFA, as opposed to y MC for prediction of MC. y FFA having passed the Fisher's F-test, further evaluation of its adequacy is performed by the residual plots, Fig. 1. Observation of the plot showed no adequate signs of block effect and the variance does not increase in trend. This is an indication that the Box-Behnken design is well randomised. Therefore, it can be inferred that y FFA can be used to predict FFA for various values of the variables within their specified minimum and maximum limits.
The resulting response surface plot for y FFA within specified limits of the process variable is shown in Fig. 2. Although y MC has not passed the F-test, its plot is also shown. Analysis of Eq. (16) and Fig. 2a, shows, CPO FFA increases with x 1 , x 2 , and x 3 , and Eq. (17) and Fig. 2b indicate MC in CPO decreases with x 2 , and x 3 .
Based on the discussion thus far, it can be inferred that in developing the CPO FFA kinetics, the effect of moisture content, and temperature should be incorporated into the model. And as would be expected temperature would be used for CPO MC kinetics. (12) k m = k m0 e ϕ�MC i /273R  Fig. 3, and Table 3, which highlights the estimated rate constant, k , evaluation criteria R 2 and RMSE for the model. Based on the R 2 value, the best curve fit was for 0.30% MC i and 35 ℃ and 0.25% MC i and 85 ℃. While for RMSE, 0.25%, 0.30% MC i and 35 ℃ was the best. In general, from the tread in the values of k , as storage time progresses, CPO FFA increases with higher MC i and temperature. This observation is also clearly indicated in Fig. 4. Figure 4, was deduced from the curve-fit of the rate constants, k , in Table 3 to Eq. (11), and the resulting Arrhenius constants,k 0 , activation energy, E , and evaluation criteria, R 2 and RMSE are given in Table 4. Analysis of the R 2 and RMSE values, indicate the model has a very good fit with the rate constants, k.
Furthermore, to account for the effect of moisture content, the deduced Arrhenius constants, k 0 , in Table 4 was curve fitted to Eq. (12) via the ratio, k m = k 0 /k 0,0.20% and MC i . The resulting Arrhenius constants, k m0 , activation energy to MC i , ϕ , and evaluation criteria, R 2 and RMSE are given in Table 5. Observation of Fig. 5, and as supported by the R 2 and RMSE values, the model shows a very good fit with k m .
The final curve-fitted model for CPO FFA kinetics is given by Eq. (18) as deduced from Eq. (10) and (13), noting that the point of initialisation for the model is 100%.    19 for FFA kinetics of extracted oil from stored almonds. However, in Lin et al. 19 report, first-order reaction kinetics were proposed with consideration of the influence of relative humidity (RH), and temperature. The reported data indicated that FFA formation was more temperature-dependent and the reaction rate increased faster at higher RH than at lower RH. The RH in Lin et al. 19 report can be assumed to have the same effect as the initial moisture content, MC i used in this work, since RH was kept constant in this case. The inference on the effect of temperature and RH in Lin et al. 19 report can be observed from the increasing rate constant, k with temperature in Table 3 and k 0 with MC i in Table 4 respectively. Furthermore, it should be noted that as opposed to the consideration of modelling rate constant as a function of initial moisture content, i.e. k 0 = f (MC i ) and keeping the activation energy, E constant in the Arrhenius model, Eq. (11). Lin et al. 19 estimated both variables as a function of RH.
Reaction kinetic for CPO MC. Similar to the CPO FFA kinetics, CPO MC kinetics were also investigated. The resulting curve-fit of experimental data to the proposed kinetic model, Eq. (10), showed a first-order reaction kinetic, n = 1.0 fits the model best. The comparison of experimental data and kinetic model are shown in Fig. 6,      Fig. 6. The reason for this poor result can be attributed to the fact that during the drying of moisture in liquid samples, the moisture is constantly in equilibrium with the moisture in the surrounding air. Therefore, if the rate of moisture removal is not high enough (i.e., at low temperature) small amount of moisture is removed, while at a higher temperature a lot of moisture is removed. And in between low and high temperatures, the rate of moisture removal is quite inconsistent. In addition, factors such as amount and effective or exposed surface area of samples, and variation of the relative humidity of the surrounding can significantly affect the rate of moisture removal. Furthermore, the rate constants, k , in Table 6 was curve-fitted as a function of temperature using the linearised form of Eq. (11), and the results are given in Table 7 and based on the R 2 value the model showed poor fit, as also illustrated in Fig. 7. In general, the poor fit of the CPO MC kinetic model, Eq. (10) and (11), collaborate with the result of lack-of-fit of the statistical multi-regression model, Eq. (17) earlier developed.
The curve-fitted CPO MC kinetics from the combination of Eqs. (10) and (11) is given by Eq. (19) with an initialisation of 100%.
By comparison of the CPO MC kinetics with the Thin-layer drying curve equations, Eq. (19), when algebraically integrated, is equivalent to the Henderson and Pabis model 36 . Having earlier highlighted the inadequacy of CPO MC kinetics, there is a need to deduce an elaborate model to predict the changes of CPO MC with temperature, following comprehensive approaches described in literature [37][38][39] .

Conclusion
In conclusion, the result of this work showed mathematically how moisture content, x 1 , temperature, x 2 , and storage time, x 3 influenced changes in CPO FFA from the analysis of a statistically significant multi-regression model, deduced from the Box-Behnken design of experiment for these process factors. It was found that each of these factors significantly influenced CPO FFA linearly ( x 1 , x 2 and x 3 ), quadratically ( x 2 1 , x 2 2 and x 2 3 ) and by an interaction between moisture content and temperature ( x 1 x 2 ) only. Additionally, a well-fitted half-order reaction kinetic model was developed to describe the changes of CPO FFA with the incorporation of the effect of temperature and moisture content via the Arrhenius model. Similarly, a multi-regression and first-order kinetic model was developed for CPO MC. However, the multi-regression was statistically insignificant, and the firstorder kinetic model had a poor fit with experimental data. The inadequacy of CPO MC models developed may suggest that elaborate considerations of physical and thermodynamic phenomena associated with the evapouration of moisture from CPO needs to be further investigated.

Data availability
The authors do not wish to share supplementary data, because all vital experimental data are the same as the data used for the plots in this article.