Optimization of phenol adsorption from synthetic wastewater by synthesized BiFeO3 perovskite material, using split-plot central composite design

Optimization of adsorption processes using statistical methods of experiment for the removal of pollutants from wastewater, in an effort to curb the global problem of water pollution, is increasingly being adopted because it is cost-effective and time-saving. In most cases, standard central composite designs (CCDs) are often employed for the optimization processes, where the experimental variables are often randomized completely. However, most experiments especially within the industries often involve factors with some hard-to-change (HTC) levels and some with easy-to-change (ETC) levels, in which case the HTC factor cannot be completely randomized, and this challenge can only be overcome by the use of a split-plot CCD. However, there is scarcity of literature on the use of split-plot CCD for the optimization of adsorption processes, and hence in this study, the prime conditions for the removal of phenol onto BiFeO3 from synthetic wastewater were studied. The effect of three adsorption variables (pH, adsorbent dosage, and shaking time) was investigated using split-plot CCD. pH was considered as the HTC factor due to the amount of time, acid and/or base required to change it, while the adsorbent dosage and contact time were the ETC factors. Quadratic model was developed for the phenol percentage removal. The optimum adsorption conditions obtained from the study were adsorbent dosage of 0.60 g, pH of 7 as well as contact time of 167 min with desirability of 1. The predicted and experimental values obtained were 89.73 and 89.21%, showing good agreement between the experimental value and those predicted by the quadratic model. Langmuir isotherm model was found to be the best fit for the equilibrium adsorption data giving rise to monolayer adsorption capacity of 106.50 mg/g. The pseudo-second-order kinetic model’s correlation coefficient (R2) was higher than that of the pseudo-second-order kinetic suggesting the applicability of the model to the adsorption of phenol. The synthesized BiFeO3 could be considered as a viable alternative to the expensive commercial activated carbon for the removal of phenols in wastewater, and the use of split-plot CCD model makes the experiment much easier to run and save time and/or cost due to fewer number of runs and restriction in the randomization of HTC factors.

and population increase which have resulted in the generation and discharge of organic and inorganic pollutants into the environment leading to the contamination of most water supplies and posing a threat to the environment, human health, and ecological system (Awad et al. (2020)). The environmental effect of organic pollutants has been recognized for decades, according to Awad et al. (2020). Organic pollutants are mostly produced by the chemical, pharmaceutical, and petroleum processing sectors, with phenol and phenolic wastes being the most common pollutants (Ruifang et al. (2020) and Garba et al. (2016)).
According to Medjor et al. (2015), phenol has the chemical formula C 6 H 5 OH and is an aromatic, weak acidic organic molecule. It is a flammable white crystalline substance, and there are about 2000 varieties of phenolic compounds in nature which can be found in sewage, natural water, and drinking water. Sabrina et al. (2019) stated that phenol contaminant has high toxicity even at low concentrations and is considered one of the most hazardous organic pollutants. Phenol inhalation, ingestion, or skin adsorption can cause coma, convulsion, cyanosis and affect the liver, kidney, lung, and vascular system. However, according to Abdelwahab and Amin (2013) phenolic compounds are hazardous and mutagenic at high concentrations. Further, Muhammad et al. (2019) stated that ammonia and phenol are among the chemicals that cause cancer, respiratory disorders, fibrosis, and pneumonia. Also, Karunarathne and Amarasinghe (2013) stated that phenol overexposure affects the central nervous system; in such situations, the results are catastrophic because the damage is irreversible, so the US Environmental Protection Agency (USEPA) has limited the permissible discharge level of phenol to less than 1 mg/L in order to minimize the hazard caused by phenol, while the limit in drinking water is 0.002 mg/L (Sunil et al. 2013). Therefore, phenol must be removed from wastewater before disposing it into the environment.
Over the years, several methods have been used for the treatment of wastewater containing phenol and other organic pollutants. Bioremediation has been proven to be a promising and ecologically friendly treatment approach of organic pollutants in soil and water (Zhao et al. 2020). Other methods used to remediate wastewater over time, according to Karunarathne and Amarasinghe (2013), include oxidation, precipitation, ion exchange, photocatalytic degradation, solvent extraction, and adsorption.
In adsorption process, activated carbon is the most used adsorbent (Lin and Juang 2009). This material is said to have good surface properties, such as a large surface area and a large pore volume. Activated carbon may potentially contain functional groups that interact with contaminating molecules. These characteristics make it a good material for adsorption activities, as it has a high adsorption capacity (Bhatnagar et al. 2013). However, there is a lot of interest in research into alternative and low-cost precursor materials for activated carbon synthesis, and perovskites nanoparticles have recently been employed as adsorbents. Perovskite nanoparticles have recently gained the interest of researchers due to their outstanding properties and promising applications in several applications such as photocatalysis, sensors, catalysts for CO 2 conversion, dry reforming catalysts, syngas generation catalysts, and recently adsorption (Garba et al. 2019).
Response surface methodology, RSM, according to Myers et al. (2009), has continued to play vital roles in developing, optimizing, and improving processes, particularly with several input variables; it is an area of experimental design which consists of a group of mathematical and statistical techniques used in the development of an adequate functional relationship between a response of interest, y, and a number of associated control (or input) variables denoted by x 1 , x 2 x 3 …, which potentially influence some performance measure or quality characteristic, of the process under study. However, one difficulty in applying classical response surface designs is that they inherently assume that all factors are equally easy to manipulate, thereby allowing for complete randomization of experimental run order. In practice, most industrial experiments cannot be completely randomized due to the presence of some factors with hard-to-change (HTC) levels (Yakubu and Chukwu 2018). Therefore, these experiments are often conducted in a manner that restricts the randomization, which leads to a split-plot structure. In a split-plot design, experimental runs are performed in groups, where in a group, the levels of the HTC factors are not reset from run to run. This creates dependence among the runs in one group and leads to clusters of correlated errors and responses. The first paper to exclusively focus on conducting response surface experiments with a split-plot structure was by Myers and Lentner, as quoted by (Yakubu and Chukwu 2018). The authors investigated efficiency of various second-order response surface designs when conducted with a splitplot structure. This was then followed by other authors like Vining and Kowalski as quoted in Myers et al. (2009), who modified completely randomized central composite and Box-Behnken designs (CCD and BBD) to accommodate a split-plot structure. Split-plot central CCD consists of factorial (f), whole-plot axial (α), subplot axial (β), and centre (c) points. The split-plot design, when compared to the completely randomized design, would provide more precise results at a cheaper cost (Arun et al. 2017). Thus, in this study, we synthesized the perovskite BiFeO 3 using the solution combustion method and investigated its optimum phenol adsorptive potential by applying split-plot central composite design in the hope of identifying a better, less expensive, and environmentally acceptable material for organic pollution remediation. To the best of our knowledge, there was no literature describing the use of the stated perovskite for phenol removal from wastewater.

Preparation of BiFeO 3 perovskite material
BiFeO 3 perovskite nanomaterial was prepared using the solution combustion method as reported by Penalva and Lazo (2018). In comparison with other physical procedures, the solution combustion method requires less time, and it is simple to use and produces a high synthesis rate for the target compounds, demonstrating superior performance in material fabrication . The method, according to Dwita and Ismojo (2016), is more simple, energy efficient, cost-effective and requires lower process temperatures than other methods. First, 200 ml solutions of the respective precursors were prepared by weighing and dissolving the required amounts (obtained from the equation of the reaction below) of the compounds in little amount of distilled water. Few drops of nitric acid were added, and the solutions were stirred using magnetic stirrer until a homogenous mixture was observed, and then, more distilled water was added to make it up to the 200 ml mark of the volumetric flask. Then, the prepared solutions of Bi(NO 3 ) 3 .5H 2 O and Fe(NO 3 ) 3 .9H 2 O were combined at room temperature with magnetic stirring, followed by the addition of glycine and urea fuels in the appropriate stoichiometry as given in the equation of the reaction. After that, drops of dilute nitric acid were added and the mixture stirred until a clear solution was obtained: The entire mixture was placed in a porcelain dish and placed on a heating mantle and gradually raised the temperature to 400 °C. The mixture was heated until all of the solvent had evaporated, and a combustion reaction occurred with the release of substantial volume of gas. The brown powder obtained was then placed in a muffle furnace and subjected to a thermal treatment at 500 °C for 1 h with a 2 °C/min ramp. The powder obtained with thermal treatment was allowed to cool and was taken for its respective characterization and use.

Characterization of the BiFeO 3 perovskite
Using a powder diffractometer equipped with Cu Kα (λ = 1.5406 Å) as the radiation source, the crystal structure of the BiFeO 3 perovskite was determined. The XRD pattern was recorded in the 2θ ranges 5-80 °C. From the XRD pattern, the average crystallite size ( D ) of the BiFeO 3 powder was estimated using the Scherer Eq. (1) where λ is the X-ray wavelength (1.5406 Å), β is the full width at half maximum (FWHM), K is Scherer constant, and θ is the Bragg angle. The morphology and particle size were studied by scanning electron microscope (SEM).

Preparation of phenol stock solution
Analytical grade of phenol was used for the preparation of the stock solution by dissolving 1 g of phenol in distilled water and making it up to the mark in 1-L volumetric flask (1000 mg/L). It was stored in a brown-coloured glass container to prevent photooxidation and was changed with fresh one after one week if not consumed (Medjor et al. 2015). Synthetic adsorbate solutions of various concentrations were prepared from the stock solution by successive dilutions.

Design of experiment using split-plot CCD
Split-plot central composite designs (CCDs) of experiment in Design Expert Software 12 version computer programme, which is a subset of RSM, were used in this study (Garba and Rahim 2014) The split-plot CCD consists of factorial (f), whole-plot axial (α), subplot axial (β), and centre (c) points (Myers et al. 2009). The experiment was developed and executed as a three-factor central composite with a split-plot structure. Three factors were considered during the design process, and these factors are: pH, contact time, and adsorbent dosage. Each of the numeric factors was set to 5 levels: plus and minus alpha Page 4 of 14 Garba et al. Bulletin of the National Research Centre (2022) 46:186 (axial points), plus and minus 1 (factorial points), and the centre point. The experimental runs were divided into groups based on pH (the HTC factor), with all runs, having the same pH level being run at the same moment in one group. The split-plot CCD used is defined by operations that include 8 factorial runs, 4 whole-plot axial runs, and 4 subplot axial runs (placed a distance α away from the centre and allowing the design to rotate), and 9 centre runs (used to determine the data reproducibility and experimental error). For the HTC factor, there are a total of 8 groups. In all, 25 experiments (runs) were generated from the software for the three variables. The low, central, and high levels, designated as − 1, 0, and + 1, respectively, are defined in Table 1. Furthermore, this approach necessitates tests outside the experimental range in order to predict response functions outside the cubic domain (denoted as ± α). The ranges of these values are displayed in Table 1.

Adsorption studies
The adsorption study was carried out using a method similar to the method used by Arun et al. (2017) and Okibe et al. (2018). Batch adsorption experiments were performed at ambient temperature with initial adsorbate concentration of 40 mg/L in the absence of illumination (UV light) to avoid photooxidation. Adsorbent dosage, pH, and contact time were set as generated by the split-plot central composite design matrix in Table 2. All the experiments were run at ambient temperature using 100 ml of 40 mg/L of the adsorbate solution in a 250-ml Erlenmeyer flask, adjusting the adsorbent dosage, pH as well setting the contact time accordingly as generated by the experimental design. The solution-adsorbent mixtures were stirred using a laboratory shaker at 250 rpm. The pH of the solution was adjusted with HCl and NaOH solutions using a pH meter. At the end of each shaking, the samples were filtered through Whatman No. 1 filter paper to eliminate any fine particles. In all the experiments, blank measurements were taken. The equilibrium concentration of phenol was then determined using a double-beam UV-visible spectrophotometer at a wavelength of 270 nm.
The adsorbate removal percentage was determined using Eq. (2) where Co and Ce are the initial and final concentrations of the adsorbate (mg/L) in the solution.
The adsorption capacity was calculated using Eq. (3) (2)  where q e is the adsorption capacity (mg/g), Co is the initial adsorbate concentration of solution (mg/L), Ce is the equilibrium adsorbate concentration (mg/L), v is the volume of adsorbate solution used (L), and w is the weight of the adsorbent (g).

Adsorption isotherms
Adsorption isotherm models are commonly used to characterize the adsorption process because they show how much of the adsorbate is adsorbed by the adsorbent at a given temperature. Adsorption data are fitted to models with assumptions that are near to the real situation, because there are multiple known isotherms based on various assumptions and circumstances, to allow comparison of different materials (Awad et al. 2020).

Langmuir adsorption isotherm
Langmuir isotherm assumes that there is an infinite number of binding sites having the same affinity for adsorption of monolayer and there is no interaction between the adsorbed molecules. Based on the assumption, Langmuir isotherm is represented in the linearized form by Eq. (4) (Idowu 2015) where C e is the equilibrium concentration (mg/L), and q e is the amount adsorbed at equilibrium time mg/g. q max and K L are Langmuir constants which relate adsorption capacity and the energy of adsorption, respectively. q max and K L can be obtained from the slope and intercept of the graph C e q e against C e (Fig. 5). R L , which is an indicative of the isotherm shape and predicts whether an adsorption system is favourable or unfavourable, is an essential element of the Langmuir isotherm. R L , according to Upenyu et al. (2017), can be calculated using Eq. (5) where Ce is the initial adsorbate concentration (mg/L) and K L is adsorption equilibrium constant (L/mg) obtained from the intercept of the Langmuir plot. The value of R L indicates the type of the isotherm to be either: ii. Linear (chemisorption and physisorption taking place at the same rate) if R L = 1 iii. Favourable (physisorption is the predominant process) if 0 < R L < 1 or iv. Irreversible (non-reversible chemisorption process) if R L = 0.

The Freundlich model
The Freundlich isotherm assumes monolayer sorption with a heterogeneous energetic distribution of active sites, combined with interaction between the adsorbed molecules. The linear equation for the model is given by Eq. (6) where q e is similarly defined as that in the Langmuir isotherm, and C e is the equilibrium concentration. A plot of log 10 q e against log 10 C e gives a straight line graph ( Fig. 6) for dilute solution. The slope is 1 n , and the intercept is K . K and n are Freundlich constants; n indicates how favourable the adsorbent's adsorption capacity is, as well as the capacity of the adsorbent-adsorbate system, which is the magnitude of adsorption capacity, and it decreases with increasing temperature (Idowu 2015).

Temkin isotherm model
The Temkin isotherm model assumes that the adsorption heat of all molecules decreases linearly with the increase in coverage of the adsorbent surface, and that adsorption is characterized by a uniform distribution of binding energies, up to a maximum binding energy. In other words, the model is based on the assumption that heat of adsorption will not remain constant. It decreases due to interaction between sorbent and sorbate during adsorption process. The Temkin isotherm can be described by Eq. (7) where K T is the equilibrium binding constant (L/mol) corresponding to the maximum binding energy, b t is related to the adsorption heat, R is the universal gas constant (8.314 J K −1 mol −1 ), and T is the temperature (K).
Plotting q e versus lnC e from Eq. (7) results in a straight line (Fig. 7) of slope RT b T and intercept RT b T lnK T .

Kinetics Studies
Kinetic investigations are important for adsorption studies because they can predict the potential ratecontrolling step and the mechanism of adsorption reactions. Adsorption kinetics was used to measure the adsorption uptake with respect to time at constant concentration. The mechanism of the adsorption process was determined using pseudo-first-order and pseudosecond-order models. The Lagergren equation (Lagergren 1898) as quoted by Nyankson et al. (2020) is used to express the pseudo-first-order model, which posits that the rate of adsorption site occupation is proportional to the number of unoccupied sites and it is given by Eq. (8) where q e and q t are the amount of phenol adsorbed (mg/g) by BiFeO 3 at equilibrium and at any time t, respectively, K 1 (min −1 ) is the rate constant of pseudo-firstorder model, and t is the time (min). A plot of ln (q e − q t ) versus t (Fig. 8) is used to estimate the values of K 1 and q e from the slope and the intercept, respectively. The pseudo-second-order kinetic model can be expressed in the linear form as given by Eq. (9) where K 2 is the equilibrium rate constant of pseudo-second-order model (g/mg min −1 ). The values of q e and K 2 were determined from the slope and intercept of the plot of (t/q t ) versus t (Fig. 9), respectively.

X-ray diffractometer (XRD)
The detected peaks in the XRD patterns were identified as rhombohedra perovskite structure of BiFeO 3 with space group R-3 m and lattice parameters a = b = c = 3.9520 (Å).

Scanning electron microscope (SEM)
The BiFeO 3 images, at 100 μm and 20 μm magnifications, exhibit agglomerated particles of various sizes and shapes.

The independent variables and their levels
The low, central, and high levels of the variables are designated as -1, 0, and + 1, respectively. Tests outside the experimental range are denoted as ± α, and they help to predict response functions outside the cubic domain.
The independent variables, their levels, and responses generated from the model

Normal plot showing the relationship between the predicted and experimental data
The plot of predicted values versus actual values shows that the data points were evenly split by the 45-degree line.

Results of restricted maximum likelihood (REML) analysis with Kenward-Roger p values for phenol removal
Variance terms are grouped into two sections: a wholeplot section for the HTC factors and a subplot section for ETC factors. The p values (Prob > F) for the whole-plot and subplot terms, as a whole, were observed to be highly significant with values much less than 0.05 which is the generally acceptable alpha level.

Results of the estimates, standard errors and variance inflation factors (VIF) of parameters of the fitted split-plot CCD model
The VIF values are equal to one except for the second and the last two terms of the model (a 2 , B 2 , and C 2 ). The second factor had VIF value of 1.04, and the last two had VIF values of 1.02 each which indicates that it was moderately correlated with all the other factors, and thus, there was no case of multicollinearity in the data.

Results of the process optimization
At the optimum conditions (pH of 7, contact time of 167 min, adsorbent dose of 0.6 g as well as desirability of 1), the experimental (actual) % adsorbed values obtained in the laboratory for phenol were only slightly lower than those predicted showing good agreement between the experimental values and those predicted from the model.

3D response surfaces for interaction between the adsorption variables
Adsorbed value about 90% was achieved at pH 10 and contact time of 243 min (keeping adsorbent dosage at 0.6 g). Adsorbed value above 90% was achieved at pH 10 and adsorbent dosage of 0.4 g (keeping time at 165 min).
Removal of above 80% was achieved with both contact time and adsorbent dosage playing important roles.

Langmuir, Freundlich and Temkin isotherm plots for the adsorption process
The Langmuir isotherm model appeared to be significantly more applicable than the Freundlich and Temkin models with the highest R 2 value ( Table 6). The fitness of the Langmuir model to the adsorption process connotes that the adsorbates ion molecules from bulk solution were adsorbed on specific monolayer which is homogeneous in nature.

Plots for pseudo-first-order and pseudo-second-order kinetic models for the adsorption process
The pseudo-second-order kinetic model correlation coefficient was R 2 = 0.998 and that of the pseudo-first-order model was R 2 0.7952 (Table 8), suggesting the applicability of the pseudo-second-order kinetic model to the adsorption of phenol.

Results of the kinetic models' constants and their correlation coefficients
The q e cal. value for phenol (7.58 mg/g) for the pseudosecond-order kinetic reaction was found to be in better agreement with the q e exp. values (6.14 mg/g) confirming the fitting of the data to the pseudo-second-order model.

Characterization of the adsorbent X-ray diffractometer (XRD) analysis
The XRD pattern was recorded in the 2θ ranges 5-80 °C. Figure 1 shows the X-ray pattern of the synthesized bismuth ferrite powder. The detected peaks in the XRD patterns were identified as rhombohedra perovskite structure of BiFeO 3 with space group R-3 m and lattice parameters a = b = c = 3.9520 (Å) which corresponds to ICDD reference No. 01-072-2112. This is identical to what Omid et al. (2015) prepared. However, the commonly observed by-products like Bi 2 Fe 4 O 9 and Bi 24 Fe 2 O 39 located at 2θ = 27.97° (Yongming et al. 2011) were not detected; hence, XRD studies proved the synthesized BiFeO 3 was highly crystalline with high purity. The average crystallite size ( D ) of the BiFeO 3 powder was estimated using Scherer equation (Eq. (1)). The average crystallite size as determined by the XRD pattern was 140.20 nm. However, because the width of diffraction peaks is inversely proportional to crystal size, the detected BiFeO 3 peaks were quite sharp, showing a high degree of crystallinity of the powder and hence the huge crystal size (Dovydas et al. 2020). BiFeO 3 particle size can also change dramatically with heat treatment/temperature, according to the literature (Yongming et al. 2011).

Scanning electron microscope (SEM) analysis
The morphological examination of the synthesized BiFeO 3 was carried out using a scanning electron microscope at various magnifications (100 and 20 μm) as shown in Fig. 2 a and b for proper view of the morphology. The images exhibit agglomerated particles of various sizes and shapes caused by the production of a considerable number of gases (NO, NO 2 , CO, NH 3 , and H 2 O) during the combustion process between fuels and nitrates (Pattnaik et al. 2018). Peñalva and Lazo (2018) and Pattnaik et al. (2018) all synthesized and reported BiFeO 3 images that were similar in morphology to ours. Table 2 shows the design matrix incorporating both the HTC and the ETC variables, their ranges, and the responses which are percentage removal of phenol (Y P ).

Design matrix for the split-plot CCD, factors, and corresponding responses
Using pH (Factor a) as the HTC factor due to the amount of time, acid and/or base required to change it, contact time (B), and adsorbent dosage (C) as the ETC factors. The pH (HTC) factor was used to group the experimental runs. The grouping of HTC factor makes the experiment much easier to run. The whole-plot column has been added to keep track of HTC factor combinations, which form the "whole plots". All of the runs in any given group maintain a constant level of the HTC factor. As can be seen from the results, the percentage removal of phenol using BiFeO 3 as adsorbent is between 66.55 and 89.21%. The application of split-plot RSM resulted in the following regression Eq. (10) for phenol, which explains the functional relationship between removal percentage and process variables where (Y P ) is % of phenol adsorbed by BiFeO 3 , pH (a), time (B), and adsorbent dosage (C), respectively. Factors in capital letters indicate ETC factors, while the small letter represents the HTC factor. The goodness of fit of the model was verified using determination coefficient (R 2 ) where the present model explained over 92% of the variation in the observations (R 2 = 0.9248) for phenol. Figure 3 shows the plot of predicted values versus actual values, respectively. The plot helps us to detect observations that were not well predicted by the model and in this case, it shows that the data points were evenly split by the 45-degree line. (10) Page 8 of 14 Garba et al. Bulletin of the National Research Centre (2022) 46:186

Restricted maximum likelihood (REML) analysis with Kenward-Roger P values for phenol removal by BiFeO 3
According to Arun et al. (2017), to get the proper P values from a split plot, specialized statistical tools such as restricted maximum likelihood (REML) must be performed. The purpose of maximum likelihood estimation is to find the parameter values that make the observed data most likely. This is another way to estimate variances. In the split-plot case, REML estimates the group variance for the whole-plot factors and the residual variance for the subplot factors. Once the variances are estimated, generalized least squares (GLS) is used to estimate the factor effects and Kenward-Roger's method is then used to produce F-tests and the corresponding p values.
Results of restricted maximum likelihood analysis for the model are given in Table 3 for phenol. We observed from this table that variance terms are grouped into two sections: a whole-plot section for the HTC factors and a subplot section for ETC factors. This result is consistent with the findings of Yakubu et al. (2020) in their experiment. The p value (Prob > F) for the whole plots for phenol was observed to be highly significant with values much less than 0.05 which is the generally acceptable alpha level. That is, for all the terms making up the whole-plot (HTC) portion of the model, it was also observed that the p values for the individual whole-plot term (a) were significant with value much smaller than 0.05. This indicates that the pH (factor a) has significant effect on the adsorption of the adsorbate. Also, the quadratic effect of this factor (a 2 ) on the adsorption of phenol is highly significant. Table 3 also shows the subplot terms as a whole were observed to be highly significant, though some individual subplot terms (aC, BC, and C 2 ) were not. The main effect of factors B and C (contact time and adsorbent dosage); the interaction effect of pH and time (aB) and pH and adsorbent dosage (aC); and the quadratic effect of time (B 2 ) were each observed to be highly significant on the adsorption of phenol (with p values far less than 0.05). Also, the computed regression coefficients and variance inflation factors (VIF) as given in Table 4 were considered. The VIF reflect the extent to which multicollinearity is present in a regression analysis by measuring how much the model's variance is inflated by a lack of orthogonality in the design (Yakubu et al. (2020)). Similar results were obtained by Yakubu et al.    (2020) in their experiment on optimization of cake height to achieve desired texture. The VIF values in Table 4 are equal to one except for the second and the last two terms of the model (a 2 , B 2 , and C 2 ) which indicates the factors had non-correlated orthogonal relationship with all the other factors in the model. The second factor had VIF value of 1.04, and the last two had VIF values of 1.02 each which indicates that it was moderately correlated with all the other factors (Yakubu et al. (2020)). Thus, there was no case of multicollinearity in the data. These statistics signify that the fitted model was good; that is, the model captures most of the variation.

Process optimization
Optimum conditions for adsorption are very vital for batch adsorption processes, and therefore, the optimization of phenol adsorption onto the synthesized BiFeO 3 was carried out by using design-expert software (version 12). In the optimization analysis, the target criterion was set as maximum values for the response (Garba et al 2016). The optimum adsorption conditions obtained were pH of 7, contact time of 167 min, adsorbent dose of   Page 10 of 14 Garba et al. Bulletin of the National Research Centre (2022) 46:186 0.6 g as well as desirability of 1 as presented in Table 5. The model validation results indicate that at the optimum conditions the experimental (actual) % adsorbed values obtained in the laboratory for phenol were only slightly lower than those predicted showing good agreement between the experimental values and those predicted from the model. The use of the model will save time because fewer experimental runs will be required in the laboratory (Okibe et al 2018).

Interactive relationships between independent variables and the response
To assess the interactive relationship between independent variables and the response, 3D surface response were utilized. As shown in Fig. 4a, % adsorbed value above 90% can be achieved at pH 10 and contact time 243 min (holding adsorbent dosage at 0.6 g). It can be observed that the improvement of removal efficiencies for phenol is attributed to a simultaneous increase in both contact time and pH. The removal of phenol as a function of pH and the adsorbent dosage is given in Fig. 4b where % adsorbed value above 92% can be achieved at pH 10 and adsorbent dosage of 0.4 g/100 ml (keeping time at 165 min). From the 3D surface, it can be noticed that the improvement of removal efficiency for phenol is attributed to the increase in pH of the mixture. It is observed that the phenol removal was not significantly affected with increasing the adsorbent dosage from 0.37 to 0.83 g/100 ml. Although the increase in surface area and adsorption sites of the adsorbent can be attributed to increasing the adsorbent dose, the values of phenol removal did not vary considerably as the adsorbent dosage was increased. The major reason could be that during the adsorption reaction, the adsorption sites remain unsaturated, whereas increasing the adsorbent dosage increases the number of sites available for adsorption. The agglomeration of the adsorbent can also be a reason for the lack of increase in phenol removal capacity at the high concentrations of the adsorbent (Salem et al. 2014). This means that the higher values of phenol adsorption are obtained by increase in the pH. The optimum adoption at high pH may be because at low pH, most of the phenol ions are in protonated state and the adsorbent surface tends to be partially positively charged as suggested by Buhani et al. (2018). This results in the occurrence of electrostatic repulsion so that the adsorption is weak. Hence, at higher pH, there is a contribution of excessive OH − free ions which are attracted to the adsorbent and tend to make it negatively charged causing the electrostatic interaction between the phenol and the adsorbent. Figure 4c shows the effect of contact time and adsorbent dosage (holding pH at 7) on phenol removal. It can be observed that removal of phenol, above 90%, can be achieved with both contact time and adsorbent dosage playing important roles (Fig. 5).

Adsorption isotherms results
The applicability of the Langmuir, Temkin, and Freundlich isotherms under the optimum conditions (as shown in Table 5) was tested using data from equilibrium studies carried out at ambient temperature with varying initial concentrations of the adsorbate. A plot of 1/q e versus 1/Ce of the Langmuir model (Fig. 6) gave a straight line with 1/q max as intercept and. 1/ K L q max as slope and hence q max and K L were estimated (Table 6). For the Langmuir model, the constants K L and q max are the adsorption equilibrium constant (mg/L) and monolayer adsorption capacity of the adsorbent (mg/g), respectively. Table 6 shows the values of these constants, which were calculated from the slope and intercept of the Langmuir plot. The q max values obtained for phenol from the linear Langmuir plot (Fig. 5) were 106.50 mg/g. The R L , which is an indicative of the isotherm shape and predicts whether an adsorption system is favourable or unfavourable, is an essential element of the Langmuir isotherm. R L , according to Upenyu et al (2017), can be calculated using Eq. 5 The Freundlich plot of log q e versus log Ce of the varied adsorption data for the adsorbate onto the BiFeO 3 adsorbent is shown in Fig. 6. The numerical values of the constants n and K f for the plot were determined from the slope and intercept of the plot and are presented in Table 6. The constant K f can be defined as the adsorption or distribution coefficient and indicates the relative adsorption capacity of the adsorbent in relation to the bonding energy. It denotes the amount of adsorbate that has been adsorbed onto the adsorbent for a given unit equilibrium concentration. For the Temkin model, plotting q e versuslnC e , Eq. (7) resulted in a straight line of slope RT b T and intercept RT b T ln K T .
The Temkin isotherm plot is indicated in Fig. 7, and its evaluated isotherm parameters are shown in Table 6. The observed large variation in R 2 from the other isotherms studied clearly showed that the Temkin equation is not well fitted with the equilibrium experimental data.
The Langmuir isotherm model appeared to be significantly more applicable than the Freundlich and Temkin models with the highest R 2 value according to the fitting findings presented in Table 6. The fitness of the Langmuir model to the adsorption process connotes that the adsorbates ion molecules from bulk solution were adsorbed on specific monolayer which is homogeneous in nature. As can also be seen from Table 6, the R L value lies between 0 and 1 which confirms the adsorption processes to be favourable (physisorption) under the studied conditions. The n value obtained from the Freundlich plot (Fig. 6) for phenol was 1.27 (greater than unity) further indicating favourable adsorption conditions as well as a physical adsorption processes for the phenol. Table 7 shows the comparison of the maximum monolayer adsorption of phenols onto various adsorbents. From the results, it can be seen that the synthesized BiFeO 3 performed better (as indicated by the q max values) than several adsorbents that were reported in the literature.

Adsorption kinetics
Adsorption kinetics measures the adsorption uptake with respect to time at constant concentration. The mechanism of the adsorption process was determined using pseudo-first-order and pseudo-second-order models. The Lagergren equation (Nyankson et al. 2020) is used to express the pseudo-first-order model, which posits that the rate of adsorption site occupation is proportional to the number of unoccupied sites as given by Eq. (8). A plot of versus t (as shown in Fig. 8) was used to estimate the values of k 1 and q e from the slope and the intercept, respectively (Table 8). On the other hand, the pseudosecond-order kinetic model was expressed in the linear Fig. 4 a 3D response surface for interaction between time and pH for phenol adsorption. b 3D response surface for interaction between pH and dosage for phenol adsorption. c 3D response surface for interaction between time and dosage for phenol adsorption Garba et al. Bulletin of the National Research Centre (2022) 46:186 form by Eq. (9), where K 2 is the equilibrium rate constant of pseudo-second-order model (g/mg min −1 ). The values of q e and K 2 (Table 8) were determined from the slope and intercept of the plot of (t/q t ) versus t, respectively (Fig. 9). The pseudo-second-order kinetic model correlation coefficient was R 2 = 0.998 and R 2 0.7952 for the pseudo-first-order model, suggesting the applicability of the pseudo-second-order kinetic model to the adsorption of phenol. The q e cal. value for phenol (7.58 mg/g) for the pseudo-second-order kinetic reaction was found to be in better agreement with the q e exp. values (6.14 mg/g) further confirming the fitting of the data to the pseudosecond-order model. K 2 value of phenol was 0.035, indicating faster adsorption rate of the adsorbate by the BiFeO 3.

Conclusions
Three adsorption parameters (pH, adsorbent dosage, and contact time) were optimized with the aid of splitplot CCD for the removal of phenol onto BiFeO 3 with the percentage removal (Y P ) as the analysis response. Based on the data analysis obtained, the individual and combined effects of the parameters were significant in the removal of the phenol adsorbate. It was established that the optimum adsorption conditions for the removal of phenol in synthetic wastewater by the adsorbent were pH 7, time 167 min, and 0.6 g adsorbent dose. The BiFeO 3 removed up to 89.21% phenol. Therefore, the synthesized BiFeO 3 could be considered as a viable alternative to the expensive commercial activated carbon for the removal of phenols in wastewaters.