1. Introduction
Wastewater treatment concurrently with biogas production via the anaerobic digestion (AD) process has been a universally adaptable technology [
1]. However, environmental pollution, water scarcity, food, and energy insecurity have become pressing matters for sustainable development in the twenty-first century [
2,
3]. Also, greenhouse gas emissions (especially CO
2) being associated with global warming, and fossil fuel combustion, undermine a sustainable environment [
4]. Herein, wastewater treatment is envisioned as a renewable energy source that can be used for biogas production as alternative sources of energy [
3,
5]. Thus, reclaiming wastewater for reuse and biogas production (biogas can be purified and used as an automotive fuel) can ease poor country’s standard of living as far as water and energy are concerned [
2].
Notwithstanding, treating industrial effluent has become extremely relevant as emerging contaminants (pharmaceuticals, pesticides, nanomaterials, etc.) and toxic chemicals originating from households and industries are posing great threats to human health and the ecosystem when discharged into the environment without proper treatment [
6,
7]. So, mitigating environmental pollution and its consequences warrant degradation of the high organic constituents of industrial effluents.
Conventionally, treatment of water and wastewater involves the use of physical, chemical, and biological methods [
7,
8]. The biological treatment involves bacterial and fungal biosorption, anoxic and anaerobic/aerobic processes [
2,
7]. Membrane filtration, coagulation-flocculation, flotation, precipitation, ion exchange, adsorption, ultrasonic mineralization, ion-pair extraction, and electrolysis [
6,
9,
10] are some of the physico-chemical treatments. Moreover, these technologies generate by-products during treatment that come with additional treatment costs [
9,
10]. Thus, it becomes very critical to design and select cost-effective technologies as an option in the water and wastewater treatment settings.
Currently, advanced oxidation processes (AOPs) are being developed to degrade nonbiodegradable contaminants [
11]. Heterogeneous photocatalysis is a unique AOP that uses photocatalysts like TiO
2 and UV light to break down pollutants into inoffensive end products like CO
2, H
2, and mineral acids [
11,
12]. The availability, low cost, and photochemical stability of TiO
2 makes it a preferred photocatalytic degradation catalyst for organic contaminants [
12]. To stimulate oxidation-reduction potential, photocatalytic reactions use photons with energy larger than the bandgap of a semiconductor, usually TiO
2 [
13]. Several researchers have investigated the photocatalytic treatment of hazardous substances in untreated wastewater [
7,
14,
15,
16]. Baseline pollutant concentration, photocatalyst concentration, and pH are among the key parameters that can influence photocatalytic activity [
13]. Other parameters including irradiation time, light intensity, light wavelength, catalyst type and temperature can also affect the degradation routes [
17]. Similarly, AD processes are influenced by several factors which includes carbon/nitrogen (C/N) ratio, temperature (mesophilic and thermophilic), organic loading rate (OLR), hydraulic retention time (HRT), and pH [
1]. This makes it difficult to assess the relative importance of more than a few influencing variables, especially when they have a multifunctional effect on the outputs.
To assess the influence of operating parameters on photocatalytic process efficiency, current research has relied on a conventional one-factor-at-time (OFAT) experimental technique [
18,
19]. While OFAT techniques are labor-intensive, they do not depict multiple variable interaction effects. Conversely, fewer runs are needed when using response surface methodology (RSM) [
20]. Box Behnken (BB), Central Composite Design (CCD), Central Composite Face centered (CCF), and full factorial are some of the most frequent design matrices used in the RSM technique. [
18,
19]. Environmental remediation studies have shown the great potential of using BBD and CCF for experimental design and optimization [
17,
18,
19,
20]. This is due to their ability to obtain more information from just a few numbers of their experimental matrix. Generally, using RSM in process optimization involves the following steps: (a) identification of response/s; (b) screening of multivariates according to the design of experiment; (c) building of an empirical response surface model; (d) and the application of various response optimizations through mathematical modeling [
17,
20].
In essence, AD (kinetically slow) and AOP (energy-intensive) processes are very complex, with limited knowledge of their integrated system (AD-AOP) [
21]. Therefore, developing the biophotocatalytic (BP) system as an alternative technology to the AD process comes in handy. Herein, the CO
2 methanation mechanism was carried out by incorporating Fe-TiO
2 into the BP system and was optimized via RSM to maximize its efficacy.
Figure 1 presents the light-driven methanation reaction scheme with the dynamic state of water-splitting and CO
2 reduction steps. Ideally, the reactant and its intermediate products emulate natural photosynthesis [
22]. Similarly, the presence of the photocatalysts being energy-driven for the initial excitation of its electrons via light absorption also promotes hydro generation potential [
21]. However, improper catalyst loading and operating conditions can limit microbial activity and biogas production, yielding poor methane quality. Furthermore, the development of the BP system is still underway, which warrants process modeling, optimization, and control to maximize the process efficiency, as well as support the lab-scale design of a pilot plant. Herein, this study employed experimental data obtained from a modified RSM- BBD matrix with input factors at three levels (−1, 0, +1) with two center points at a constant temperature of 37.5 °C and an organic loading rate of 2.38 kgCOD/Ld. This was aimed at investigating the relationship between the input variable (catalyst dosage and HRT) and the design outputs (biogas, COD, color, and turbidity) via the modified RSM-BBD. Additionally, the analysis of variance (ANOVA) was employed statistically to ascertain the significance of the response models, whereas numerical optimization was used to optimize and maximize the system desirability for degradation of the organics (COD) for the biogas production.
3. Results and Discussion
In principle, photoexcited catalysts can reduce CO
2 to CH
4 or split water molecules (H
+; O
2−) via light-induced redox reactions [
22,
26]. As a result (
Figure 1), the H
2 produced during the exergonic metabolism is subsequently utilized by the methanogens to enhance the CO
2 reduction reaction into methane [
1,
5,
27,
28]. Results of the characterized municipality wastewater showed a high organic load (2380 ± 14 mgCOD/L) with a VS/TS ratio of 0.75 (
Table 1). The VS/TS ratio > 0.5 proves the wastewater used was biodegradable [
11,
22].
Figure 3 shows the weekly monitoring of the BP system operated under anaerobic conditions, whereby the degraded organic content (2380 > 115 > 87 > 30 > 12 mgCOD/L) increased biogas production. The cumulative biogas recorded for week 1 to week 4, were 750, 1950, 1980, and 100 mL/d, respectively as depicted in
Figure 3. Also, the activated photons of the Fe-TiO
2 catalyzed the microbes to enhance the degradability of the wastewater organic content.
Figure 4 shows the increased methane level of the biogas produced recorded for week 1 to week 4, respectively, like 79%, 83%, 80%, and 95%, with the rest as CO
2 composition. From our preliminary studies, there was the need to investigate the effect of the HRT and catalyst load on the BP system efficiency. Therefore, understanding their individual or interaction effects on the BP system led to the RSM studies at a constant temperature of 37.5 °C, and an organic loading rate of 2.38 kgCOD/Ld.
3.1. RSM Modelling and Statistical Analysis
According to the modified RSM-BBD matrix (
Table 2), 13 randomized runs were carried out, including one center point (run 10) and four duplicates (runs 1:3; 2:9; 5:8; and 7:12) as depicted in
Table 3 with additional information in
Table S1. All the responses of the center point (run 10) showed little variation, indicating the consistency of the experimental runs. In all the experimental conditions, the removal efficiency of COD, color, and turbidity was observed to range from 90–98%. This degradation efficacy corresponded to biogas production of 125–335 mL/d. However, the biogas experimental results (125–335 mL/d) were estimated to deviate by 5% to that of the modified RSM-BBD model predicted results (130–327 mL/d). Subsequently, the significance of the factors was determined by modeling the experimental data as a function of the individual response (
Table 3). This resulted in a reduced quadratic model with the Design Expert software that provided equations (3–10) in their coded and actual factors for the distinctive responses. The linear (A, B), interaction (AB), and quadratic (A
2, B
2) terms represent the hierarchal form of the model, which combined the forward and backward regression option [
29]. This regression option augments and eliminates the insignificant variables that do not meet or fall short of the model levels required (
p < 0.05) [
19,
29]. The positive and negative coefficients in the predicted models (3–10) denote energetic and antagonistic effects of the process variables, respectively. In essence, the positive coefficient suggests that the response will be favored by an increase in such variable interaction. Conversely, the negative signs implies the response efficiency will decrease with an increase in those variables [
29].
3.2. Analysis of Variance (ANOVA)
The response data was analyzed with the Design expert software, and the derived models were then fitted using the ANOVA (
Table 4). The ANOVA shows how well the quadratic models fit the experimental values, parameters including F-value, probability >F, and adequate precision, which is a measure of error, or the signal-to-noise ratio were used [
19,
29]. The summary of the ANOVA for the various quadratic models is depicted in
Table 4. The models’ variables depict the two-factor interactions between the catalyst load (A) and the HRT (B) for the biogas (Y
1), COD (Y
2), Color (Y
3), and Turbidity (Y
4) efficiency. All the variables and their interactions were significant in the proposed models except the interaction of AB of the COD (Y
2) quadratic model (5–6) as indicated by a probability value being more than 0.05 (
Table 4). A significant interaction between AB means that the effect of each variable depends on the value of the other variable [
20]. Thus, increasing the catalyst load (AD) will require a longer period to reduce its potency. Adequate precision for all the models was greater than 4, suggesting a high signal/noise ratio; hence, the models can be used to navigate the design space [
29].
Subsequently, the extent of correlation was estimated using experimental-predicted data interactive plots that is represented by plotting the predicted values against the experimental ones (
Figure 5).
Figure 5a,c,d shows a strong linear correlation between the experimental and the model predicted data with their high regression coefficients (R
2) as depicted in
Table 4. However, in
Figure 5b, only a few data points were frequently scattered around the diagonal line, which could account for the on-and-off degradation of the organics (%COD removal) based on the subjected experimental conditions. The standard error (SE) for the straight line of best fit had insignificant deviation (
p < 0.05) at a 95% confidence level. As can be inferred in
Table 4, the COD (Y
2) model regression is very low, even though all the model’s Adjusted R
2 and Predicted R
2 values are in reasonable agreement with a difference of less than 0.2. The standard error (SE) for the straight line of best fit shows insignificant deviation (<0.05) at 95% confidence levels.
3.3. One-Factor-At-Time Technique
3.3.1. Effect of Catalyst Load
The impact of the Fe-TiO
2 load over the range of 2–6 g was used to study the degradability of the organic content of the wastewater into biogas. In
Figure 6, it is observed that the degradation efficiency as a function of (a) biogas, (b) COD, (c) color, and (d) turbidity increased with an increase in the catalyst load. However, after the optimum regime (>4.5 g) there was a drop in the degradation efficiency. The decrease in efficiency might be due to excess hydroxyl radicals generated which agglomerated and reduced the active surface area [
30]. In essence, increasing the catalyst load increased the collision frequency between the hydroxyl radicals and the organic content interspecies, which diminished the photoactivity of the system [
27,
30]. Also, the exponential growth of the microbial organism and their degradability activity increased biogas production [
27]. This suggests that the induction of the catalyst facilitated the conductive electron interspecies transfers, which increased the microbial activity and the methanogenic activity [
13,
30,
31].
3.3.2. Effect of Hydraulic Retention Time (HRT)
The HRT investigated from 1–31 days had a significant effect on the digester performance as it facilitated microbial activity to the organic interspecies contact time, which increased its digestibility.
Figure 7 shows an initial increase in degradation (COD, color, and turbidity) during the first 5 to 15 days followed by a decline (15–31 days), whereas a fairly constant rate for biogas production (
Figure 7a) was observed. This indicates that a longer HRT is required for the methanogens to digest the organic component and convert it to biogas [
32]. In addition, the digestion process will slow down if volatile fatty acids accumulate due to acidic bacteria being dominant with less organic content [
26,
30,
33].
3.4. Response Surface Interaction Plots
The modified RSM-BBD was used to illustrate the interactive impact of the factors on the response.
Figure 8 shows the presence of interactions between the factors, catalyst load, and HRT (AB) on (a) biogas (b) COD, (c) color, and (d) turbidity. The graphical representation of the three-dimensional (3D) surface plots of the response models was selected based on the influential factors and their interaction that can be utilized to maximize the system desirability. Curvatures of a significant magnitude can be seen in the graphs (
Figure 8). These curves indicate that the correlation between the factors (AB) and the response (biogas, COD, color, and turbidity) was well fitted on a quadratic function (3–10). This elucidates
Figure 8 as the degradation efficiency increased to the maximum with an increase in HRT; likewise, with the catalyst load (
Figure 6). However, the curvature (
Figure 8) highlights an increase in catalyst with dosage from 2 to 4 g and thereafter decreased with a further increase in catalyst load (5–6 g). The trend (
Figure 8), observed at higher dosages of the catalyst, can be attributed to particle re-stabilization [
30,
33]. In addition, an overdose and excess of hydroxyl radicals can also aggregate to decrease the active surface area [
30]. In essence, at the stage of re-stabilization, there is a charge reverse between the interspecies [
26,
30,
33]. As illustrated (
Figure 8) with an arc-line, the optimum regime can be observed within the HRT of 15–25 days. This validates the positive sign of the response quadratic models (3–10) as reported earlier, which had significant impact on the system predictability (
Figure 5).
3.5. Response Optimization and Confirmation Test
The numerical optimization technique was employed to maximize the responses (biogas, COD, color, and turbidity) and determine the optimal conditions (with respect to the experimental runs) using their respective quadratic equations (3–10). Also, the comparative evaluation of the experimental and model predicted responses at each level of the experimental conditions (
Table 3) showed considerable correlation with the two input variables (A and B) as a constraint. The goal for the optimization was defined as a function of the input variables set within the range of their levels such as low (−1) and high (+1) to maximize the desirability of the responses. The optimal solution selected out of 78 solutions (
Table S1) to maximize the responses is presented in
Figure 9 with 100% desirability. As it can be inferred from the ramp plot (
Figure 9), 267 mL/d of biogas, 97.75% COD, 98% color, and 99% turbidity removal were attained at 4 g catalyst load and HRT of 21 d. The selected optimal conditions validated and confirmed experimentally as presented in
Table 5 were in good agreement with the predicted values. This suggests the model’s predictability was consistent (
p < 0.05) at 95% confidence levels with a small standard error (SE) and standard deviation (STD) (
Table 5). Also, the results obtained were compared with other studies (
Table 6), which showed the appreciable efficacy of the BP system for wastewater treatment and biogas production.