3.1. Total phosphorus removal and hydraulic loading rate
Table 1 shows descriptive statistics for the fourteen water quality parameters monitored in the MSLs influent. For instance, the pH measured in the raw wastewater has minimum, maximum, and average values of 7.20, 8.90, and 8.17, respectively. In the same context, to examine the effect of the hydraulic shock load on MSL efficiency in removing TP, the five MSL units were subjected to an increasing HLR (between 250 and 4000 L/m2/day).
Table 1
Summary statistics of the HLR and the MSLs influent characteristics
Parameter | Unit | Range | Mean | Standard deviation |
HLR | L/m2/day | 250–4000 | - | - |
pH | - | 7.2–8.9 | 8.17 | 0.35 |
DO | mg/L | 0.36–1.3 | 0.78 | 0.31 |
EC | µS/cm | 1722–2421 | 2015 | 198 |
TSS | mg/L | 206.4–380.2 | 278.9 | 39.7 |
BOD5 | mg/L | 244–396 | 314.4 | 46.5 |
COD | mg/L | 430–619 | 504.2 | 47.1 |
NH4+ | mg/L | 15.4–34.15 | 23.1 | 4.8 |
NO2− | mg/L | 0.11–0.93 | 0.43 | 0.24 |
NO3− | mg/L | 10.3–34 | 20.5 | 5.44 |
TKN | mg/L | 22.3–49.3 | 33 | 6.3 |
TN | mg/L | 57–109.3 | 77 | 12.42 |
PO43− | mg/L | 3.1–7.2 | 4.9 | 1.3 |
TP | mg/L | 5.2–9.9 | 7.4 | 1.14 |
TC | log units | 5.9–7.4 | 6.83 | 0.47 |
Regarding the MSLs performance, we note that there are significant (p < 0.05) differences between the efficiency of the MSL systems subject to different HLRs and becomes insignificant (p > 0.05) when the HLR rises from 2000 to 4000 L/m2/day (Fig. 3a). On the other hand, obtained results show a decrease in the TP level from 7.40 mg L− 1 to 0.73 mg L− 1 for the HLR of 250 L/m2/day, and 3.96 mg L− 1 for the HLR of 4000 L/m2/day (Fig. 3a). Thus, it can be noticed that increasing the HLR value negatively affects the level of TP in the MSL effluent. This finding appears to be confirmed in the scatterplot (Fig. 3b), where we can see a linearly decreasing trend of TP removal as a function of HLR. In addition, for a better understanding of the performance and stability of the MSL systems, Fig. 3c depicts the changed curve of TP content with time increased. As illustrated, the effluent concentrations of the MSL systems (MSL1, MSL2, and MSL3) were relatively stable over the experiment duration. In contrast, the level of this indicator was unstable in the effluents of the MSL systems (MSL4 and MSL5) exposed to a high HLR (2000 and 4000 L/m2/day), which indicates that by increasing the HLR, the performance of the MSL system decreases and becomes unstable. These findings are in line with what has already been published in the scientific literature [9, 35], which found that increased HLRs may reduce the MSL system's TP removal efficiency.
Overall, despite the decreased efficiencies demonstrated by the MSL4 and MSL5, all the MSL systems consistently maintained significant (p < 0.05) TP reduction, with removal rates ranging from 47 to 90%. In addition, the HLR is found to be inversely related to the hydraulic shock load applied to the MSL systems' surface area. As the HLR increases, the wastewater's residence time within the MSL decreases, reducing their interaction with the MSL layers and, therefore, decreasing the MSL's purification efficiency [36]. This finding was also supported by Toreu et al. [37], who found that P adsorption by soils rises as reaction time increases.
Regarding P removal mechanisms in the MSL system, many studies [38–40] have reported that its removal is done by adsorption on iron hydroxides and aluminum hydroxides present in soil mixture and soil colloids. An et al. [8] have reported that iron metal in the soil mixture immobilizes P by adsorption and precipitation during the percolation of wastewater. Indeed, soil mixture iron is transformed into ferrous iron (Fe2+) and then transported to the gravel layer (aerobic conditions) to be oxidized to ferric ion (Fe3+), allowing for P purification in the MSL system once again [41]. Therefore, increasing P contact time with adsorbents inside the MSL system is a key component in its removal. However, short-cut flow or preferential flow would lead to the reduction of contact time between wastewater and the MSL media, which will negatively influence the transformation of iron to ferric ions. Furthermore, Song et al. [15] reported that MSL parameters such as aeration, pH, and wastewater dispersion can also influence TP removal. Indeed, sufficient oxygen could promote the oxidation of iron to ferric ions in the aerobic zone (gravel layers) resulting in increased P adsorption by soil particles, while pH could influence the proportions of hydroxyl ions (OH−) [42]. Jianbo et al. [43] have previously verified this assertion, stating that the pH of a phosphorus solution has a significant impact on the whole adsorption process, notably on the adsorption capacity.
3.2. Effect of climatic variables
In this paper, the importance of rainfall, temperature, and evaporation on the performance of MSL is examined. Figure 4a shows that the rainiest months are from January to April, and the driest months are from June to August. Regarding precipitation, it is heavier in the spring and winter, whereas the summer is drier with rainfall from 3 to 13 mm. Regarding temperature, the coldest month was January (2.6°C), while the hottest months are June to August (45°C), which explains the significant increase in evaporation during this period (between 207 mm and 326 mm).
In the same context, Fig. 4b depicts the efficiency of TP removal in the MSL systems throughout the experiment period. This removal was stable in the MSL1 during the period (Jul.17 - Jun.18), with an average removal between 86% and 92%. In comparison to the MSL3, which varied somewhat between a decrease (Sep.17 - Oct.17) and an increase (Nov.17 -Jun.18), the MSL2 also showed good performance for TP removal during the experiment period (between 77% and 89%). However, the performance of the MSL systems (MSL4 and MSL5), fluctuated significantly between an increase and a drop, wherein the MSL5 recorded the lowest removal in October (16.7%). This result can be explained by the low concentration of TP in the MSL influent in October (as indicated in Fig. 3). Therefore, this might imply that when TP concentrations in the influent are low, the MSL system's performance subjected to high HLR is limited.
On the other hand, Pearson's coefficient was also used to assess the relationship between TP removal and climatic variables (Fig. 4c). The results reveal that the linear correlation was weak (-0.1 ≤ r ≤ 0.03), with no significant (p > 0.05) relationship between local climatic data and the MSL performance. This result was consistent with Masunaga et al. [39] findings, which indicated that seasonal variables did not appear to have a significant impact on TP level in the MSL system. Overall, the MSL system's performance has remained constant throughout the year, with no significant influence of climatic variables in the semi-arid study area, especially for longer HLRs. However, the increase in the HLR makes its performance unstable.
3.3. Feature selection
Regarding an output variable, when the number of candidate input variables becomes large, selecting those who have significant relationships with the output could be an effective strategy [44, 45]. In this study, the Pearson coefficient and BIC metrics were used for this purpose. The pairwise correlation matrix (Fig. 5a) indicates the input variables that change linearly with TP removal. Each square includes both the correlation coefficient (r-Pearson) and the significance asterisks. We can note a significant (p < 0.05, r = − 0.77) negative relationship between TP removal and the HLR (Fig. 5a). However, this relationship is significantly (p < 0.05) positive with the level of TSS, pH, NO3−, NO2−, TP, and TC in the MSL influent, where Pearson's coefficient value ranges between 0.15 and 0.26.
Regarding the BIC metric, ten subset selections were evaluated and compared based on these candidate significant variables. This task was processed using the lmSubsets R-package (version 3.5.2). Figure 5b shows that the first subset has the lowest BIC value, while Fig. 5c gives more details about the ten subset selections. The candidate variables are indicated on the x-axis, while the y-axis shows the selection number. The results of these selections were very similar, with only a few variables switching between being added and excluded. Therefore, the BIC-selection consisting of the features (HLR, pH, PO43−, and TP) is suggested to be significantly (p < 0.05) pertinent to explain TP removal in the MSL system.
In line with this result, Lamzouri et al. [17] also used pH, TP, and HLR as input features to predict the TP content at the outlet of the MSL system. Kotti et al. [46] have also used HLR as an important input variable to predict the removal of P in free-water surface CWs. Thus, the empirical equation developed by Akratos et al. [47] was based on HLR as an important factor for predicting the removal of P in horizontal CWs. Vidal et al. [48] reported that pH appeared to be a good estimator for the reduction of TP in a sand filter. Overall, the selected variables are also used in previous studies as informative predictors related to the removal of P in wastewater treatment systems.
3.4. Tuning parameters
In this study, the goal of tuning the data-driven models is to identify a set of parameters that would result in an optimal model with lower error values and higher prediction accuracy. Thus, to tune these solutions, the parameters of each model were assessed, and the optimal combination was selected based on the highest R2 value. Figure 6 illustrates the tuning parameter profile for each method. In the case of the NN model (Fig. 6a), accuracy increases in parallel with the increase in the number of hidden neurons (size = 6) and reaches the highest R2 value (R2 = 0.90) for six neurons for weight decay of 0.03. Thus, the appropriate NN model architecture to predict TP removal in the MSL system is 4 units in the input layer, 6 units in the hidden layer, and one unit in the output layer (Fig. 6b). Regarding the RF tuning profile (Fig. 6c), as the number of mtry increases, the accuracy increases considerably (R2 = 0.93), showing that four mtry appear to achieve better performance for the RF model. Furthermore, determining the number of trees that makes the error rate stable during the tuning process is an important step.
As can be seen in Fig. 6d, the model error decreases as the number of trees rises. Therefore, the optimal number of trees to use, according to RF's final model, is 500. Regarding the KNN model, as was already noted, the grid search approach was employed to determine the final values of the Minkowski distance and the maximum number of neighbors (kmax) while maintaining the kernel function as "optimal" for weighting. Thus, the accuracy increases in lockstep with distance (Fig. 6e). Therefore, the set of parameters (kmax = 3, distance = 3, kernel = optimal) provides the best performance for the KNN model.
In the same context, after training the models on 80% of the total data set, they were compared and evaluated using three performance metrics. A comparison of the estimated performance is shown in Fig. 7a. In the case of the MLR, the model shows modest performance, followed by the KNN model. In terms of the MAE, its mean value varies from 9.2% (MLR) to 4.64% (RF). Furthermore, the mean RMSE values are roughly consistent with the MAE values: 11.96% (MLR), 8.34% (KNN), 6.55% (RF), and 7.26% (NN).
Overall, the RF model was the most accurate data-driven model during the training process. These conclusions are supported by the boxplots in Fig. 7b, which show that the residual distributions for all the models are different. Furthermore, the plot (Fig. 7b) suggests that the residuals of the RF model are close to those of the NN model, but are typically less than those of the KNN and MLR models.
3.5. Data-driven model evaluation
After tuning the models and identifying which one performs well, the performance metrics were used to assess the data-driven model's prediction accuracy. Figure 8 illustrates a comparison between the predicted and the real values for the TP removal in the MSL system using the validation data set (20%).
Regarding the MLR model (Fig. 8a1), some real data points are overestimated, resulting in a significant difference between real and predicted values. This might be owing to the MLR model's failure to account for the non-linear relationship between the predictors and the output variable. In addition, the model’s ability to predict TP removal is unsatisfactory in that the RMSE, MAE and R2 are respectively 12.12%, 9.22% and 0.64 (Fig. 8a2). However, although the KNN model presents an improvement over the linear model, some data points are underestimated (Fig. 8b1). This result can be confirmed by the fit between the real and predicted values in that the RMSE, MAE and R2 were 7.18%, 4.61% and 0.87 respectively (Fig. 8b2).
In the case of the RF, it can be assumed that this model shows a nearly perfect match between the real and predicted values (Fig. 8c1). Furthermore, the predictive performance of the RF model (RMSE = 5.29%, MAE = 3.63% and R2 = 0.93) in this article outperformed the previous models in terms of TP removal predictions (Fig. 8c2). Similarly, for the NN model, Fig. 8d1 visualizes a good convergence between the predicted and real values. This finding can be confirmed by the good fit between the real and model predicted values (R² = 0.90). In addition, the estimated values of RMSE (6.11%) and MAE (4.51%) demonstrate that the NN model was also suitable for predicting TP removal in the MSL system (Fig. 8d2).
In summary, the MLR and KNN models have exhibited modest to desirable prediction accuracy, while the RF and NN models showed good prediction accuracy. These findings support the predictive power of the NN in predicting MSL system performance. However, the RF model was shown to be more accurate in this investigation. Furthermore, as compared to previous studies [15, 17] whose accuracy in predicting TP removal in the MSL system was ranged between (R2 = 0.60) and (R2 = 0.85), the RF model developed in this study improved this prediction (R2 = 0.93) based on relevant and informative input variables.
Finally, due to the number of parameters to tune (mtry) and the obtained accuracy, as well as simplicity in the implementation, it can be concluded that RF is a simpler and more powerful data-driven approach that accurately predicts TP removal in the MSL system. Furthermore, as compared to the SCA and NN methods, the balance between complexity and accuracy of the RF method makes it a decent choice for predicting MSL performance.
3.6. Feature importance
Following the prediction of TP removal, feature importance analysis was performed to explore the most important variable that influences this removal. Based on the accurate model, the importance of each predictor was calculated using variable_importance( ) function. Figure 9 indicates that all of the input variables selected using the feature selection technique are important and explain the output variable.
The results demonstrated the HLR's supremacy as the most important parameter controlling TP removal in the MSL system, as well as the significance of their nonlinear relationship. Furthermore, TP level in the MSL influent is ranked as the second most important variable, followed by PO43− and pH, respectively. Thus, these findings are consistent with the literature [9, 46, 48], implying that the level of pollutants in raw wastewater, in combination with the hydraulic shock load, governs MSL removal efficiency.