Modeling and optimization of process parametric interaction during high-rate anaerobic digestion of recycled paper mill wastewater using the response surface methodology

One of the most critical industries globally, especially in Malaysia, is the recycled paper industry. There are 20 paper mills in Malaysia, 18 of which use raw materials of 100% recycled paper. From this recycled paper industry, recycled paper mill wastewater (RPMW) is being generated. RPMW is categorized as complex wastewater due to significant amounts of suspended solids, biological oxygen demand (BOD), chemical oxygen demand (COD), lignin, ammonia, volatile fatty acid (VFA), and biodegradable organics (Zwain et al. 2016a). Therefore, RPMW must be treated according to permissible limits before being released into a municipal sewer system or receiving streams.

Owing to high amounts of readily biodegradable materials, RPMW could be anaerobically treated with different combinations of biological systems (Cai et al. 2019). Anaerobic digestion (AD) is one of the most appropriate and efficient methods that have been used to treat industrial and municipal wastewater. It has the potential to produce methane and a semi-stabilized digest, which can be further processed as a substitute for inorganic fertilizers (Chatterjee & Mazumder 2019). However, for the effective utilization of AD in industrial wastewater treatment, the successful development of high-rate anaerobic reactors is required. By maintaining the biomass in the reactor for a longer time, more microorganisms will be available in the system, resulting in a higher degradation process of wastewater. In this regard, the anaerobic baffled reactor (ABR), a modification of an up-flow anaerobic sludge blanket reactor, was found to be an effective system for various wastewater treatments. It is a staged reactor that promotes biomass retention within the system by driving water flows up and down across several compartments (up to eight) (Zhu et al. 2015). Subsequently, the efficacy of AD is defined by the microorganism behavior and quantity, which are influenced by the wastewater composition, system configuration, operation condition, and the microbial content formed in the anaerobic reactor (Zwain et al. 2019a). To further validate the efficacy of AD and the operational performance of the reactor system, statistical modeling and optimization studies are essential for long-term operation to verify its practicability in the RPMW treatment.

Therefore, the study's main objective is to model the performance of a high-rate novel ABR (HR-NABR) treating RPMW using the response surface methodology (RSM). This novel system is designed to overcome the drawbacks of the conventional ABR by modifying the reactor structure and microorganisms' growth mechanisms through (1) changing the baffles’ hydraulic surface to improve the contact between the organisms and RPMW; (2) by considering the slow growth rate of methanogenesis compared to acidogenesis bacteria, the volumes of the compartments are customized to have different hydraulic retention times (HRTs) to offer a suitable condition for the co-growth of acidogenic and methanogen bacteria; and (3) packing materials were placed in the system to combine the development of anaerobic attached and suspended growth microorganisms. Specifically, the study aims to evaluate the parametric interaction of influent chemical oxygen demand (COD_in) and HRT with 10 responses of COD and BOD removal efficiencies, COD and BOD removal rates, lignin removal, methane content and production, effluent pH, effluent alkalinity, and effluent VFA using a three-level factorial design (3²-FD). The predicted results were validated, and the optimal operating conditions of the HR-NABR for achieving a high operational performance during the RPMW treatment were then determined.

RPMW was used as a substrate in this study and collected from the Muda paper mill in Simpang Ampat, Penang, Malaysia. It had an average COD concentration of 3,812 mg/L, a BOD concentration of 1,868 mg/L, a volatile suspended solid concentration of 1,967 mg/L, and a biodegradability (BOD/COD) ratio of 0.49. The system used is an HR-NABR operated at a laboratory scale, and the schematic diagram is shown in Figure 1. It has five compartments separated by modified vertical baffles with a total effective operational volume of 35 L. The system is designed to combine the growth of suspended and attached microorganisms by adding packing materials to the system. Further information about the system configuration was previously illustrated (Zwain et al. 2017). The reactor was inoculated with flocculant anaerobic sludge supplied from an anaerobic pond treatment system treating the palm oil mill effluent from Malpom Palm Industries Bhd, Penang, Malaysia.

The reactor start-up process was previously carried out (Zwain et al. 2016a, 2016b). After that, the system process performance of the HR-NABR was investigated by varying COD_in concentrations (1,000–4,000 mg/L) and HRT (3–1 day(s)) in nine transition steps, 2 weeks each, for a total period of 126 days (Zwain et al. 2018), as shown in Table 1. Different COD concentrations were obtained by diluting the RPMW. At a constant HRT, the reactor was first operated at COD_in of 1,000 mg/L, then increased to 2,500 and finally 4,000 mg/L. Subsequently, the reactor HRT shifted from 3 to 2 days and finally 1 day, varying COD_in from 1,000 to 2,500, then 4,000 mg/L. At the end of each operational experiment, COD and BOD removal efficiencies, COD and BOD removal rates, lignin removal, methane content and production, effluent pH, effluent alkalinity, and effluent VFA were analyzed.

The RSM was used to model the HR-NABR performance and design a series of 13 experiments, and experimental design conditions are presented in Table 2. The 3²-FD was selected to analyze the interaction of two independent variables (COD_in and HRT) on the HR-NABR performance responses. Ten responses of COD and BOD removal efficiencies and rates, lignin removal efficiency, methane content and production, effluent pH, effluent alkalinity, and effluent VFA were evaluated. The chosen independent variables (i.e., COD_in and HRT) were studied at three levels, i.e., low (−1), central (0), and high (+1), and the parameter ranges are shown in Table 3. Of 13 experiments generated, nine were organized in a factorial design. The other four were to assess the replication of the central point to get a reasonable estimation of the experimental error.

The accuracy of predicted values over experimental values was evaluated using RSM analytical tools. Model fitness was tested by assessing the importance of determination coefficients of R², predicted R², and adjusted R². As much as R² is close to 1, the response reaches its ideal value and falls within the perfect intervals. Additionally, the predicted R² must reasonably agree with the adjusted R² by not more than 0.2. Adequate precision measures the signal-to-noise ratio. A ratio of 4 and more is desirable to indicate acceptable model discrimination. The model has a strong enough signal to be used for optimization.

To obtain an efficient HR-NABR performance, an optimization of experimental responses was performed using the Design Expert Software. As independent variables, high COD_in concentration and a short HRT were selected as a goal for parametric optimization. In contrast, maximum COD, BOD, and lignin removal efficiencies, high methane production, and low-effluent VFA were set as a goal response. The optimum condition was numerically generated, and it was determined with the highest overall desirability.

Predicted data over experimental responses are shown in Table 4, and the model's desirability coefficients are shown in Table 5. Good statistical models of best fit have been established to have R² values between 0.75 and 1 (Ghaleb et al. 2020). Regression data showed that the coefficient of determination R² is of close fit to 1 for most of the responses and slightly low for BOD removal (0.78), effluent alkalinity (0.82), and methane content (0.84). The high R² value indicates that the model can better estimate responses (Gopal et al. 2021). A lower close fit of predicted data could be due to low variation in responses to a change in operational variables resulting from the nature of the biological process reaching a steady state. Additionally, the values of adjusted R² and predicted R² were within 0.2, except for effluent alkalinity. However, adequate precision measures the range of predicted response relative to its associate error (Jalaludin et al. 2016). The results proved that an adequate precision value was >4 for all responses; the lowest was 7.9 for effluent alkalinity, yet it confirms that this model is desirable. As suggested by Torabi Merajin et al. (2014), low values of the coefficient of variation (0.77–22.59) showed good reliability and precision of the experiment's raw data.

The ANOVA analytical results for all responses are shown in Table 6. As different responses were studied, various polynomial models were used to fit the experimental data (responses). To evaluate the curvature effects, the experimental data (responses) were fitted to a higher degree of polynomial equations (i.e., quadratic, modified quadratic, two-factor interaction (2FI), and linear). In the Design Expert software, the experimental data (responses) were processed by default. Some experimental data may not be fitted, and transformation, which uses a mathematical function to some response data, might be necessary to meet the assumption that makes the ANOVA accurate. As errors (residuals) were a function of the response's quantity and to obtain the best data fitting, inverse square root and natural log transformations were required for the effluent alkalinity and methane production responses, respectively. The transformation type was selected based on the lambda value suggested by the Box–Cox plot. The chosen model terms are obtained after the elimination of insignificant variables and their interactions. In this study, a confidence level of 95% was selected to analyze the experimental results. P-values and F-values obtained the significance of each model. The P-values of <0.05 reveal that the model terms are statistically significant. The model modification was carried out by reducing insignificant terms with P>0.05 maintaining the hierarchy order.

Generally, the smaller the P-values and the larger the magnitude of F-values, the higher the significance of the correlating model. In this study, the P-values of the models were 0.0007, 0.0005, <0.0001, <0.0001, <0.0001, <0.0001, 0.0045, <0.0001, 0.0007, and <0.0001 for COD removal (%), BOD removal (%), COD removal rate (g/L day), BOD removal rate (g/L day), lignin removal (%), effluent pH, effluent alkalinity (mg/L), effluent VFA (mg/L), CH₄ (%), and CH₄ production (L CH₄/d), respectively. According to the statistical results, all models were significant with low probability values ranging from 0.0045 to <0.0001. It is observed that the models of independent variables were significant at 95% confidence level. The models’ -F-values were 15.7; 17.48; 183.81; 204.46; 34.84; 1,607.81; 9.13; 84.76; 15.56 and 62.49 for COD removal (%), BOD removal (%), COD removal rate (g/L day), BOD removal rate (g/L day), lignin removal (%), effluent pH, effluent alkalinity (mg/L), effluent VFA (mg/L), CH₄ (%), and CH₄ production (L CH₄/d), respectively. The model F-value of >4.5 implies that the models are significant. The interaction term of AB was insignificant for the equations defining COD and BOD removal rates (g/L day), effluent pH, effluent alkalinity (mg/L), and effluent VFA (mg/L). Detailed results of the response models are shown in the following sections.

The P-values obtained from the ANOVA in Table 6 for COD_in and HRT were 0.348 and 0.0002, respectively. These values indicated that the individual effect of COD_in was non-significant, but the impact of the HRT was highly significant. This was also supported by the fitted correlations from Equation (2), whereby a comparison between coefficients of COD_in and HRT (i.e., 0.006 and 17.3, respectively) indicates that a long HRT has higher favorite effects on COD removal efficiency than a decrease in COD_in. However, Equation (2) shows a non-significant interaction between COD_in and HRT (AB). For BOD removal, the individual P-values from the ANOVA analysis were 0.0006 and 0.0083 for COD_in and HRT, respectively. These values indicated that the effects of both factors were significant. As shown in Equation (3), the relevant coefficients of COD_in and HRT were 0.003 and 1, respectively, which reveals that an increase in COD_in and HRT have almost similar effects on BOD removal efficiency. Equation (3) also shows that there is no interaction effect between COD_in and HRT on BOD removal. This might be due to the change in BOD removal efficiency being in the range of 93–99%, whereby the model cannot explain the 6% variation in BOD removal efficiency.

The P-values of <0.0001 for COD_in and HRT reveal significant individual effects on the COD and BOD removal rates. From Equation (4), the coefficients of COD_in and HRT parameters in COD removal were 0.003 and −1.14, respectively. Also, according to fitted correlations from Equation (5), the coefficients of COD_in and HRT were 0.0006 and −0.85, respectively. These indicate that an increase in COD_in has similar effects to a decrease in HRT. Nevertheless, their interactions (AB), P-value (<0.0001), and coefficients of 0.00025 (Equation (4)) and 0.00015 (Equation (6)) confirmed the significant factorial interaction between COD_in and HRT. From the ANOVA results in Table 6, the individual effects of COD_in and HRT on lignin removal were significant model terms evidenced by P-values of <0.0001 and 0.0007, respectively. Also, the coefficients of 0.004 and 11 (Equation (6)) for COD_in and HRT, respectively, indicated that the COD_in had a higher impact on lignin removal than the HRT. However, Equation (6) articulates that the COD_in and HRT interaction is statistically non-significant.

A quadratic model was chosen to explain the response surface of pH level to changes in COD_in and HRT. The P-values obtained from the ANOVA analysis in Table 6 were <0.0001 for both COD_in and HRT, indicating that the single effects of COD_in and HRT were enormously significant. This was also supported by the coefficient of determination shown in Equation (7). There was also a strong interaction between COD_in and HRT, which was also noticed by the corresponding P-values (<0.0001) and a correlation coefficient of −6.8×10⁻⁵.

Regarding the effluent alkalinity, data transformation of inverse square root and a modified quadratic model were needed for fitting. From the fitted correlations in Table 6, it can be seen that COD_in and HRT were statistically insignificant. The coefficients (Equation (8)) of both COD_in and HRT variables in effluent alkalinity correlation were −0.043 and 6.62×10⁻⁶, respectively. This means an opposite effect for COD_in and HRT. However, the interaction of COD_in and HRT was significant because the correlation coefficient is −4.6×10⁻⁶. In addition, the relevant coefficient of P-value <0.0231 in the correlation of effluent alkalinity confirms the significant interaction between COD_in and HRT. From the ANOVA results in Table 6, the individual effects of COD_in and HRT on effluent VFA were significant with a P-value of <0.0001 for both variables. Furthermore, their interaction (AB) was also significant with a P-value of 0.0056. In magnitude viewpoint, the adjacent effects of COD_in and HRT were presented by the relevant coefficients of 0.01 and −3.3, respectively. The coefficient of interaction between them is −0.002, which expresses the minor significance of interactions of these two parameters.

As shown in Table 6, the modified quadratic model reveals that the COD_in has a significant effect (P-value=0.0135) on methane content, while the impact of HRT was insignificant (P-value=0.3983). The interaction between COD_in and HRT was found to have insignificant effects on the CH₄ content. From Table 6, the impact of COD_in and HRT was significant with a P-value of <0.0001. Nevertheless, the correlation coefficients (natural log) of 0.0005 and −1.12 for COD_in and HRT (Equation (11)) showed slightly significant model terms. However, their interaction (COD_in and HRT) was statistically non-significant.

The interaction of COD_in concentrations and HRT on the HR-NABR performance responses was analyzed using the 3²-FD. Experimental results of the COD_in and HRT interaction on 10 responses of 13 designed experiments are shown in Table 4.

Figure 2(a) shows the concurrent effect of the COD_in and HRT on the COD removal efficiency achieved from Equation (2). At a long HRT of 3 days, increasing the COD_in from 1,000 to 4,000 mg/L has increased the COD removal efficiency from 93 to 96%. Under a HRT of 1 day, the COD removal efficiency slightly decreased when the COD_in was increased from 1,000 to 4,000 mg/L, but this could be due to the effect of short HRT, not COD_in increment. At constant COD_in of 4,000 mg/L, as the HRT was shortened from 3 days to 1 day, the COD removal dropped from 96 to 83%. These results indicate that the HRT had more impacts on the COD removal efficiency compared to the COD_in. An increase in the COD removal efficiency is due to the efficient consumption of available substrates by microorganisms. In contrast, a change in parametric conditions may result in temporary microorganisms' shock that requires some time for acclimation to the new loading condition to achieve stable performance (Zwain et al. 2018). The effects of COD_in and HRT variables on the BOD removal efficiency are shown in Figure 2(b). The BOD removal efficiency was stable at the range of 93–99% and slightly affected by any change in COD_in and HRTs. This revealed that the HR-NABR performance was very stable in handling biodegradable substances.

According to standard B Malaysian Department of Environmental (Department of Environment (DOE) 2009), the discharged RPMW should contain <250 mg/L of COD and 50 mg/L of BOD. In this study, the effluent COD concentration did not exceed the limits and ranged from 64 to175 mg/L, except for experiments with an OLR >2 g/L day that have COD_in=4,000 mg/L and HRT=2 days, COD_in=2,500 mg/L and HRT=1 day, and COD_in=4,000 mg/L and HRT=1 day, whereas the COD effluent concentrations were 529, 412 and 690 mg/L, respectively. Nevertheless, an effluent BOD of 59 mg/L exceeded the limits at only the highest OLR of 4 g/L day (COD_in=4,000 mg/L and HRT=1 day) and was within the limit range of 14–48 mg/L for the rest of the experiments. Figure 2(c) and (d) shows the effects of COD_in and HRT on COD and BOD removal rates. Moreover, the removal rates of COD and BOD increased with an increase in the COD_in concentration and a decrease in the HRT. At COD_in of 4,000 mg/L and HRT of 1 day (OLR of 4 g/L day), about 3.265 and 1.875 g/L day were reached at a steady state for COD and BOD removal rates, respectively. The 3D surface plot curvatures reveal that the effect of COD_in was slightly less than a short HRT.

Figure 2(e) illustrates the change of the lignin removal efficiency as a response to the two factors investigated. Increasing the COD_in from 1,000 to 4,000 mg/L at any constant HRT has improved the lignin removal efficiency. Likewise, the lignin removal efficiency increased, at any COD_in, when the HRT was further shortened from 3 to 2 days. This means that COD_in positively impacted the lignin removal, while HRT has a favorite effect at a long HRT of 3 and 2 days only. At the HRT of 2 days, the highest lignin removal was achieved due to balanced lignin feed and microbial community development. In this multimicrobial community HR-NABR system (Zwain et al. 2017), the AD of organic compounds is accompanied by the slow degradation of lignin. However, the lignin removal decreased when the HRT shortened further to 1 day, and the COD_in increased to 4,000 mg/L. In general, RPMW contains complex lignocellulosic substances, where its lignin content degradation is slow under a high rate of AD. At a short retention time, native lignin is reluctant to enzymatic attacks by anaerobic microorganisms; however, lignin degradation compounds can be slowly degraded in complex microbial communities (Khan & Ahring 2019).

In the case of effluent pH, alkalinity, and VFA, a strong correlation was observed between the responses, not only with the parametric variation. Figure 3(a) shows the effects of the COD_in and HRT on effluent pH; as the variable increased, two opposite impacts of the COD_in and HRT on effluent pH were observed. At a long HRT of 3 and 2 days, a change in COD_in slightly affected the effluent pH level and remained >6.5. In contrast, the HRT had a reverse effect on the effluent pH; the lowest effluent pH of 6 was observed when the HRT shortened to 1 day but slowly recovered when the COD_in increased from COD_in 1,000 to 4,000 mg/L. The pH level is the critical factor during the continuous operation of an anaerobic digester. In an anaerobic system, the optimum pH condition for anaerobic microbial activity is neutral pH. Reduction in the pH level is due to the fermentation of organics to organic acids by acidogenesis. A stable pH level >6.5 is due to the efficient consumption of VFAs (Zwain et al. 2019b).

The effect of COD_in and HRT on effluent alkalinity is presented in Figure 3(b), and Table 4 shows that different values were obtained when variables changed. At a long HRT of 3 days, an increase in the COD_in had a favorable effect on effluent alkalinity. Likewise, shortening the HRT to 2 days significantly improved the effluent alkalinity. A higher HRT resulted in additional buffering that increased the metabolic activity in the system, leading to the efficient conversion of organic acids by methanogenesis (Zwain et al. 2018). Conversely, a reverse impact on effluent alkalinity was obtained at a shorter HRT of 1 day and when the COD_in increased from 1,000 to 4,000 mg/L at HRT of 2 days and 1 day, respectively. This could be due to the faster production of organic acids and their slow conversion to methane (Vuitik et al. 2019). These results indicated that both HRT and COD_in have a strong influence on alkalinity production. However, the overall effluent alkalinity proved the balance between acidogenic and methanogenic processes.

Moreover, the simultaneous effects of COD_in and HRT on effluent VFA are shown in Figure 3(c). With COD_in increasing from 1,000 to 4,000 mg/L, the effluent VFA production also increased at all HRTs, indicating a positive correlation of the COD_in on effluent VFA. This might be due to the high activity of fast-growing acidogenesis compared to methanogenesis. In addition, the concentration of effluent VFA was slightly lower at a HRT of 3 days; it was roughly increased at shorter HRT conditions, revealing that HRT had affected the effluent VFA positively. Low-effluent VFA concentration at a longer HRT can be caused by the low OLR and balanced acidogenesis and methanogenesis. Also, a short HRT supported the accumulation of intermediary products (VFA) (Leitão et al. 2006). Specifically, the three responses were strongly correlated, whereas an increase in the VFA effluent concentration was associated with a lower effluent alkalinity and pH. Nevertheless, the accumulation of VFA at a low HRT also resulted in low COD removal (83%), as shown earlier in Figure 2(a). Also, high alkalinity production was correlated with the high production of biogas in the system. When the OLR increases, fast-growing hydrolytic and acidogenic microorganisms ferment organic substrates to organic acids, leading to VFA accumulation and pH reduction (Roshanida et al. 2018).

Figure 4(a) shows the effect of COD_in and HRT on methane content. By increasing the COD_in from 1,000 to 4,000 mg/L, the methane content increased at all HRTs. This proves the favorite effect of COD_in on methane content. Regarding the HRT variable, two opposite trends for methane content were obtained at different HRTs. At constant COD_in, the methane content was increased when the HRT was shortened from 3 to 2 days, and then methane content decreased as an opposite trend when the HRT was further reduced to 1 day. A gradual increment in the OLR may enhance the system's activity by allowing a sufficient substrate for microorganisms to digest, new biomass generation, and methane production. In contrast, a high HRT will increase the substrate available for fast-growing acid-producing bacteria that inhabit methanogenesis (Hassan et al. 2014).

Figure 4(b) illustrates the simultaneous effects of COD_in and HRT on methane production. The methane production responded differently to the two variables. At a long HRT of 3 days, the methane production increased from 0.082 to 0.237 L CH₄/g COD when the COD_in was increased from 1,000 to 4,000 mg/L. A major increase in methane production was noticed when the HRT was decreased from 3 to 2 days at the COD_in of 1,000 mg/L. The high methane production is due to active methanogenesis's optimum development that efficiently converted VFAs to methane (Zwain et al. 2017). At the HRT of 2 days, when the COD_in was increased from 1,000 to 4,000 mg/L, the methane production decreased from 0.259 to 0.173 L CH₄/g COD. At the HRT of 1 day, the methane production slightly increased from 0.168 to 0.179 L CH₄/g COD when the COD_in was increased from 1,000 to 4,000 mg/L. The reduction in methane production at a short HRT might be due to inhibition of slow-growing methanogenic bacteria by excessive generation of VFA (Zwain et al. 2018).

These results indicate positive and reverse impacts for the COD_in and HRT at different conditions, and no specific trend was observed. The methane production might be better understood if it is explained in the OLR instead of HRT and COD_in. Furthermore, the CH₄ content and production fluctuated throughout the experimental runs, immediately dropped after each change in the COD_in and HRT, and gradually increased to its stable value at the end of each experiment. Regarding the interaction of methane production with other responses, slow-growing methanogenic bacteria require some time and system buffering capacity by alkalinity to consume high VFA accumulation and stabilize the system. This finding was also associated with an increase in alkalinity in the reactor. This signifies the high biodegradability of RPMW in the HR-NABR. Thus, methane production from anaerobic reactors depends on pH, alkalinity, and VFA stability (Chen et al. 2020).

The goal responses of maximum COD, BOD, and lignin removal efficiencies, high methane production, and low-effluent VFA were numerically analyzed using the RSM to define the optimal operational parameters of the HR-NABR. The optimum HR-NABR performance conditions could be associated with better effluent quality and sufficient balance between acidogenesis and methanogenesis indicated by high organics and lignin removal, high methane production, and low VFA discharge. With desirability up to 0.68 obtained by the Design Expert Software, the optimum condition for the HR-NABR in treating RPMW was obtained at the HRT of 2 days and COD_in of 4,000 mg/L. At the optimum condition, removal efficiencies of 91, 98, and 70% were obtained for COD, BOD, and lignin, respectively, whereas 26 mg/L of VFA and 8 L CH₄/d (corresponding to 0.125 LCH₄/g COD_removed) were produced.

The optimum conditions were close to the highest removal efficiencies of 96% for COD, 99% for BOD at a COD_in of 4,000 mg/L, and HRT of 3 days. At the same time, up to 71% of lignin removal was achieved at a COD_in of 4,000 mg/L and HRT of 2 days. Also, the highest effluent alkalinity of 785 mg CaCO₃/L was attained at a COD_in of 1,000 mg/L and HRT of 2 days (corresponding to an OLR of 0.5 g COD/L day). Furthermore, the highest methane production of 20.8 L CH₄/d (corresponding to 0.259 L CH₄/g COD) was produced at a COD_in of 1,000 and HRT of 2 days, while methane content of 85% was attained at a COD_in of 2,500 mg/L and HRT of 2 days. It is believed that, at the optimum condition of HRT of 2 days and COD_in of 4,000 mg/L, there was a partial separation between acidogenesis and methanogenesis, where the available substrate fermented to organic acid, and VFAs efficiently converted to methane. Before this point, there was no sufficient substrate to generate active microorganisms. Beyond that point, further substrate may cause an organic shock and inhibit the balance between acidogenesis and methanogenesis.

Table 7 compares the HR-NABR performance with other treatment systems treating RPMW in Malaysia (Mohamad et al. 2008; Muhamad et al. 2013; Muhamad et al. 2015; Dahlan et al. 2020). Owing to insufficient data for CH₄ production rate and lignin removal, the treatment system's performance was compared in terms of COD removal only. This study noticed that the HR-NABR has the highest COD removal efficiency among all the treatment systems with relatively high methane generation of 8 L CH₄/d (corresponding to 0.125 LCH₄/g COD_removed) was predicted at optimized conditions. Meyer & Edwards (2014) reviewed that for influent COD of 600–1,500 mg/L, anaerobic treatment of RPMW typically resulted in 58–86% COD removal, whereas methane generation could only be achieved at 0.24–0.4 L/g COD removed. In this study, as a comparison, the HR-NABR presented satisfactory COD removal efficiency for a wide range of RPMW feeding COD at a high HRT.

Operational parameters are the main factors influencing the performance of bioreactors. Interaction between these operational parameters is essential for optimizing the RPMW treatment process. In this study, the RSM proved to be a powerful tool for modeling and optimizing the performance of the HR-NABR treating RPMW. The models for all variables were significant with P-values of <0.0007. The interaction between COD_in and HRT was only significant for COD and BOD removal rates (g/L day), effluent pH, effluent alkalinity (mg/L), and effluent VFA (mg/L). Besides, the optimal conditions were established at 4,000 mg/L COD and HRT of 2 days, corresponding to predicted COD, BOD, and lignin removal efficiencies of 91, 98, and 71%, respectively, and VFA effluent concentration of 26 mg/L and 0.125 L CH₄/g COD_removed methane yield.

The authors declare there is no conflict of interest.

The authors sincerely acknowledge the RU-I grant scheme (A/C. 1001/PJKIMIA/814148) and iconic grant scheme (A/C. 1001/CKT/870023) awarded by the Universiti Sains Malaysia. The partial financial support of Al-Mustaqbal University College is highly appreciated.

All relevant data are included in the paper or its Supplementary Information.