The soil and groundwater samples used in this study were collected from an abandoned pesticide factory in Jiangsu Province, China. Benzene, toluene, ethylbenzene, and xylene (BTEX) was found to contaminate the core area of the site. The soils used in this work were taken from the non-polluted area at the edge of the site, and the uncontaminated groundwater was obtained from the upstream of the contaminated site. The soil samples were dried (for a week) in a greenhouse before being screened, and collected with a particle size of less than 2 mm for later use.

Benzene was added to the collected uncontaminated groundwater to a final concentration of 0.25 mM (equivalent to 20 mg/L). The groundwater was then aliquoted to the standard 40 mL threaded vials filled with the pre-prepared soil and capped with a Teflon/silicone gasket to prevent the volatilization of contaminants. PS (purity >98%, Sinopharm Chemical Reagent Co., Ltd., China) and nZVI (Sinopharm Chemical Reagent Co., Ltd., China) were added to the above vials according to experimentally designed concentrations of RSM-BBD and sodium hydroxide solution (0.1 M) and hydrochloric acid solution (0.1 M) were used for pH adjustment. All the microcosmic vials were kept for 12 h at the design temperature after uniform oscillation. Then the sample were filtered through a 0.45 μm membranes before analysis. The residual benzene concentration analysis was conducted using a gas chromatography-mass spectrometrometry system (GC-MS, 7890A/5975C, Agilent, United States) equipped with a capillary column (J&W Scientific DB-624 60 m × 0.25 mm × 1.4 μm, Agilent) after a purge and trap concentrator (Eclipse 4,552and4,660, OI Analytical, United States). Three parallel experiments were set in each group to eliminate experimental errors and determine the reproducibility of the results, and one blank control group was set without degradant to eliminate interference caused by adsorption. The removal rate of benzene was calculated by the following formula: Benzene removal efficiency = C 0 − C C 0 × 100 % (1) C 0 and C were the concentrations of benzene in the system at the beginning and the end of the reaction.

Design Expert software (DES, Version 11.0) was used to analyze the influence of different environmental parameters on the degradation rate of benzene (response) in contaminated groundwater. A four-variable with three-level BBD experiments with three central points were designed to investigate the degradation trend of benzene in contaminated groundwater in the PS/nZVI oxidation system. Four independent variables were selected, including 1) oxidant (PS) dose, 2) activator (nZVI) dose, 3) pH, and 4) temperature, with three levels designed for each independent variable (coded values + 1, 0, and −1, see Table 1). The BBD-RSM method was used to design 29 sets of experiments under different process conditions (Table 2). These 4 independent factors were presented as X 1 , X 2 , X 3 , and X 4 , respectively for statistical computations. A quadratic model was established to fit the experimental results and show the relationship between all involved environmental variables. The quadratic model is as follows: Y = b 0 + b 1 X 1 + b 2 X 2 + b 3 X 3 + b 4 X 4 + b 12 X 1 X 2 + b 13 X 1 X 3 + b 14 X 1 X 4 + b 23 X 2 X 3 + b 24 X 2 X 4 + b 34 X 3 X 4 + b 11 X 1 2 + b 22 X 2 2 + b 33 X 3 2 + b 44 X 4 2 + ε (2)

In the formula, Y is the response factor predicting the degradation efficiency of benzene, b 0 is the fixed response value at the design centre point, b i (i = 1, 2, 3, 4) represents the linear coefficient, b ii (ii = 1, 2, 3, 4) is the quadratic coefficient, b ij (i = 1, 2, 3, 4, j = 2, 3, 4) represents the effective coefficient of interactive regression, X i , X j represent independent variables, X ij represents the effect of interaction between independent variables, X i 2 and X j 2 represent the secondary effects of each independent variable, ε is the statistical error. Analysis of variance (ANOVA) was used to analyse the established model. The significance of each coefficient in the equation was determined by the F test and p-value, and the prediction ability of the model was evaluated through graphical analysis and numerical analysis.

The learning network used in this research was a hierarchical feedforward neural network with a back-propagation method, also known as an improved BP learning algorithm. Levenberg-Marquardt (LM) algorithm was selected for supervised learning functions to train the network because of its advantages of fast convergence and high computational accuracy (Srivastava et al., 2021a). The four variables determined by the response surface method were used as input layer neurons (04), and output layer neuron (01) was set as the response of benzene degradation rate. This research adopted the double-layer structure network with better generalization ability than the single-hidden layer structure which optimized the network structure of the hidden layer (Figure 1). Estimation of the number of hidden layers and the number of neurons in each layer was the main difficulties in optimizing the structure of the neural network, this problem was solved by applying the thumb rule (Gazzaz et al., 2012). In all sample datasets, 70% were used for training the network learning, 15% were used for forming the verification set, and the other 15% were used for forming the network test set (Taqvi et al., 2017), to cross-validate the ANN model for benzene degradation by the PS/nZVI oxidation system. The specific model parameters in the ANN were shown in Table 3. The scenario settings used to optimize the hidden layer structure of the ANN model were shown in Table 4. The development of the ANN model was completed by MATLAB 9.8.0.1323502 (R2020A).

The prediction performance of the models was evaluated via error functions that include the coefficient of determination ( R 2 ), root mean square error (RMSE), mean absolute deviation (MAD), and mean absolute percentage error (MAPE), expressed as: R 2 = 1 − ∑ I = 1 N A t − F t 2 ∑ I = 1 N A m − F t 2 (3) MAD = ∑ i = 1 n A t − F t n (4) RMSE = ∑ i = 1 n A t − F t 2 n (5) MAPE = ∑ i = 1 n A t − F t A t n × 100 (6) Where A t is the predicted value obtained by ANN or RSM, F t is the experimental/observed value, A m is the average of the predicted value, and n is the number of samples.

Under the optimum conditions given by RSM and ANN, the degradation kinetics of benzene with different concentrations (0.0625–1.25 mM) by the PS/nZVI oxidation system was studied. Degradation experiments were carried out on benzene solutions of specified concentrations based on optimized process parameters, and samples were collected at 20-min intervals to measure the residual benzene concentration in the system. The degradation results of benzene in the optimized experiment were modeled using the Pseudo-first order (PFO) and Pseudo-second order (PSO) kinetic models, which were widely used in the analysis of the degradation mechanism (Chowdhury et al., 2016). The diffusion mechanism of the degradation system was studied by the intraparticle diffusion (IPD) kinetic model (Zulfiqar et al., 2019).

i. PFO kinetic model

PFO kinetic model is an extensively adopted model to simulate the degradation process, which could be represented as follows: log q e − q = log q e − K 1 t / 2.303 (7)

In the formula, K 1 (L/min) represents the degradation rate constant, q e (mg/g) represents the degradation amount of benzene when the system reaches equilibrium, q (mg/g) represents the degradation amount of benzene at any time point in the reaction process, and t represents the time (min). The Curve Fitting Tool of MATLAB 9.8.0.1323502 (R2020A) was used to obtain the parameter values.

ii. PSO kinetic model

Five different forms of PSO kinetic models were used to fit the degradation rate results of benzene, namely,: Type-1 PSO, Type-2 PSO, Type-3 PSO, Type-4 PSO and Type-5 PSO. These Lagergren equations are extensively used for the degradation of liquid-solid phase based on solid capacity, and these PSO kinetic formulas could be expressed as following equations (Zulfiqar et al., 2018): t q = 1 K 2 q e 2 + 1 q e t  Type − 1  PSO (8) 1 q = 1 K 2 q e 2 1 t + 1 q e  Type − 2  PSO (9) 1 t = K 2 q e 2 q − K 2 q e 2 q e    Type − 3  PSO (10) q t = K 2 q e 2 − K 2 q e 2 q q e  Type − 4  PSO (11) 1 q e − q = 1 q e + K 2 t  Type − 5  PSO (12)

Where q e (mg/g) and q (mg/g) respectively represent the amount of benzene removed by the oxidant per unit weight under equilibrium and at any time in the reaction process. K 2 (g/mg min) represents the PSO chemical degradation rate constants. Using Type-1 PSO, the graph against t versus t/q was plotted. The slope and intercept were used to achieve the values of constant factors. In term of Type-2 PSO kinetic model, the slope and intercept were used to receive the values of constant factors after plotting the graph against 1/t versus 1/q. In the case of Type-3 PSO kinetic equation, the graph between 1/q versus 1/t was plotted and obtained the constant values using slope and intercept. By using a Type-4 kinetic model, the graph plotted against q/t versus q for obtaining the values of constant factors from slope and intercept. Similarly, by employing the Type-5 PSO kinetic equation, the graph was plotted against t versus 1/(q_e-q) and calculate the values of K₂ and q_e by using the slope and intercept.

iii. IPD kinetic model

Another kinetic theory considered the controlling factor of the reaction rate of the PS/nZVI oxidation system from the perspective of solid-liquid mass transfer. The degradation process of the benzene molecules in the system could be summarized into three steps: 1) the mass transfer of the pollutant through the outer boundary layer at the solid-liquid interface; 2) the diffusion of the liquid through the solid-phase particle channels; and 3) the oxidative degradation of the SO 4 · − on the surface of the activator (Ocampo-Perez et al., 2011). The mass transfer rate or particle diffusion rate, or the coupling effect of the two mechanisms affected the overall degradation rate of nZVI-activated PS. The IPD kinetic model obtained the prediction of the reaction process through the correlation between t 0.5 and q t , which was given by the following equation: q t = K id t 0.5 + C i (13)

Where K id (mg/g min 0.5 ) and C i refer to the rate constant of IPD and the boundary layer thickness of solid particles in the system, respectively.

The BBD-RSM successfully simulated and predicted the PS/nZVI oxidation of benzene in synthetic contaminated groundwater at each design level of the four factors. The benzene degradation rate was between 53.52% and 98.29% in the batch experiments (Table 2). The ANOVA of BBD-RSM for benzene degradation with the PS/nZVI technology showed that the predicted degradation rate of the quadratic model was highly consistent (p < 0.0001) with the actual degradation rate (Supplementary Table S1).

The value of determination coefficient R² (0.948) indicated that the suggested quadratic equation was useful for predicting the benzene degradation rate in the PS/nZVI oxidation system within the range of experimental conditions. The R²-adj value (0.895) of the RSM model was close enough to the R² value (0.948), indicating that the precision of the suggested response model was only slightly affected by the insignificant model terms. In general, R² values would increase by adding an item which has insignificant effects (p-value>0.05) to the proposed model, but the adjusted R² (R²-adj) representing significant effects would not increase (Salarian et al., 2016). In addition, the predicted R² (0.792) also has favorable anastomose with the R²-adj (0.895) (a difference of less than 0.2 was considered acceptable), further indicating that the constructed BBD-RSM model had acceptable reproducibility.

According to the screening principle of the significance of RSM model variables (p-value <0.05 and F value >5), independent variables X 1 , X 2 , X 3 , X 4 , and X 3 2 were believed to have the most significant contributions to the accuracy of the suggested model and played the most important role in the simulation results. X 1 had the largest F value (146.20), indicating that its influence on the model was dominant compared with other factors. The measured results were fitted using a quadratic model, to achieve the following regression equations: Y = 9.08 + 0.758 X 1 + 0.280 X 2 − 0.481 X 3 + 0.174 X 4 − 0.150 X 1 X 2 + 0.008 X 1 X 3 + 0.056 X 1 X 4 + 0.080 X 2 X 3 − 0.012 X 2 X 4 + 0.132 X 3 X 4 − 0.041 X 1 2 − 0.152 X 2 2 − 0.322 X 3 2 − 0.158 X 4 2 (14)

The prediction accuracy of the proposed model for the degradation of benzene by the PS/nZVI oxidation system was verified using four different evaluation methods. The measured value and predicted value of the degradation rate of benzene showed a high degree of consistency (Supplementary Figure S1A), which provided the most direct evidence for the great performance of the prediction model. The externally studentized residual analysis showed that the error values were normally distributed along a mathematical expectation value μ) close to zero and a constant variance ( σ 2 ), which confirmed the adequacy of the proposed model (Supplementary Figure S1B). In addition, the residuals of the experimental values and predicted values were uniformly distributed within a rectangular region centered at zero, indicating that the random error distribution of the proposed model was uniform and reasonable, which indirectly proved the stability of the RSM model (Supplementary Figure S1C,D).

The factor interaction analysis of the response surface methodology showed that the interaction between the PS and nZVI dosages exerted the dominant effects on the benzene degradation in synthetic contaminated groundwater (Figure 2A). In the response surface of the interaction between PS and nZVI dosages, the degradation rate of benzene showed a drastic increase with the dosage of PS increasing from 0.25 mM to 1.25 mM (53.52%–98.29%), while the effect of the activator on the reaction result was more moderate. However, when the amount of activator alone was controlled to increase from a 2:1 M ratio to 6:1, a 15% increase in benzene degradation rate was also observed (Figure 2A). This may owe to the SO 4 · − played a major role in benzene degradation, and the interaction between oxidants and activators directly determined the rate of the SO 4 · − generation (Lominchar et al., 2018). Whereas, the effect of the activator on benzene degradation was much weaker, probably due to the activation potential not fully released at the set ratio with the oxidizer (Xu Z. et al., 2021). In conclusion, the PS and nZVI dosages were positively correlated with the benzene degradation rate, although PS dosage showed the more crucial role.

The response surface of the interaction between the PS dosage and pH showed that the degradation rate of benzene decreased with the transition from the acidic to the alkaline solution at a constant PS dosage (Figure 2B). Although both H + and OH − have an activating effect on PS, alkaline conditions inhibited the transformation of nZVI to the activating Fe 2 + , which suppressed the activation capacity of nZVI in the system. In addition, OH − in the system would coordinate with Fe 2 + and Fe 3 + to form compounds with no activation efficiency, and reduce the activation ability (Rodriguez et al., 2014). Therefore, under acidic conditions and high PS dosage, the PS/nZVI oxidation system was more efficient for the benzene degradation. The response surface under the interaction between the PS dosage and temperature (Figure 2C), nZVI dosage and pH (Figure 2D), nZVI dosage and temperature (Figure 2E), and temperature and pH (Figure 2F) were relatively flat or behaved as the waterfall, which indicated that the interaction of these dependent variables had no significant effect, or was mainly caused by a single factor on the benzene degradation rate in synthetic contaminated groundwater.

The optimization results of the topological structure of the ANN model showed that when both hidden layers contain eight neurons, the evaluation function MSE, RMSE, MAD, and MAPE have the smallest values (Table 4, Scenario BP5), which represents a perfect match between the predicted value and the actual value. However, if the number of neurons in any hidden layer exceeds 8, the error of the ANN model increased instead. It was speculated that too many neurons in the hidden layer may lead to overfitting, increase the error of the test set and further increase the overall error function (Mutasa et al., 2020). Therefore, the optimal topology structure of the ANN model for the PS/nZVI oxidation of benzene was [4-8-8-1], in which case the lowest performance evaluation functions value and the highest R 2 value (0.9801) were obtained.

The simulation of benzene degradation rates using the ANN model structured in scenario BP5 showed high goodness of fit for training, validation, and test subset (Figures 3A–C), with the R 2 values of 0.9995, 0.9791, and 0.9725, respectively. In the overall model, the ANN model fits well with the “perfect fit line” (Figure 3D), and the R 2 value (0.9801) was slightly improved compared with that of the RSM model (0.948). Therefore, the developed cascade forward ANN model can be successfully used to simulate and predict the PS/nZVI oxidation of benzene in contaminated groundwater.

The empirical modeling tool ANN displayed better prediction ability than RSM in the simulation and prediction of the PS/nZVI oxidation process of benzene in contaminated groundwater (Figure 4). We compared the observed values with the predicted values of the quadratic model obtained by RSM and the trained ANN model, respectively. Their predicted performance parameters were analyzed and the deviation of the calculated values from the actual values of the two models was plotted. Compared with the RSM, the predicted value distribution of ANN was closer to the actual value. The evaluation function values for ANN were lower, and the values of the R 2 were closer to 1, which were signs of a well-fitted model (Supplementary Table S2). The highly predictive performance of the ANN was due to its ability to extract the basic interaction between dependent variables and independent variables with high accuracy without considering the degree of nonlinearity between variables, whereas the RSM only allowed fitting data based on mathematical equations (Lopez et al., 2017). It is generally accepted that RSM is usually used for the whole process of industrial system design starting from the experimental design. Its advantage reflects in putting forward more credible suggestions on the optimization of process parameters (Jiang et al., 2020). By contrast, ANN may be more suitable for processing massive experimental data, and focus on providing and establishing more detailed interrelationships between independent variables and dependent variables (Kasiri et al., 2008). Therefore, ANN may be a more powerful and flexible empirical modelling tool for remediation simulations of groundwater contaminated by organic compounds.

Three groups of process parameters recommended by the proposed models that should completely degrade benzene in synthetic contaminated groundwater were selected and the predictive ability of the proposed models was externally validated by replicated experiments (Supplementary Table S3). The results showed that under the condition of complete degradation of benzene predicted by the RSM model, the measured degradation rates of the three experimental groups were 100%, 99.87%, and 98.79%, respectively, and the predicted benzene degradation rates by the trained ANN (BP5, which with the best predictive performance) were 99.86%, 99.98%, and 99.10%, respectively. The agreement between the experimental and predicted values indicated that it was feasible to use the proposed model to simulate the actual degradation situation. The MSE between the predicted and measured values of the two models were 0.4937 (RSM) and 0.0426 (ANN), respectively, which indicated that the simulation errors of both models were acceptable in the external validation phase, but ANN still outperformed RSM.

The kinetic parameters ( q e , K 1 and K 2 ) and correlation coefficients ( R 2 ) of benzene degradation in synthetic contaminated groundwater by the PS/nZVI oxidation system were determined by linear regression method (Table 5). The R 2 values of different types of kinetic models were distributed between 0.3982 and 0.9995, and all types of kinetic models fit the experimental data well except for very few scenarios (e.g., Type-1 and Type-4 PSO at 10 mg/L benzene). The linear fitting diagram of the kinetic model showed that the fitting curve of the PFO deviated significantly from the experimental results (Figure 5A), and the parameter values of the PFO model also showed no significant regularity during the increase of benzene concentration from 0.0625 mM to 1.25 mM (Table 5), indicating that the PFO kinetic model may not be a reasonable model to explain the benzene degradation process. Type-2 and Type-3 PSO kinetic models showed the best fitting with the observed data of benzene degradation within the focused benzene concentration range of 0.0625–1.25 mM (Figures 5B,C), with R 2 values exceeding 0.99. The relationship between benzene degradation concentration versus time was more appropriately explained using these two types of PSO kinetic models, and they were fitted using 1/t versus 1/q or 1/q versus 1/t as the independent and dependent variables, respectively. Meanwhile, the instantaneous reaction rate models for different benzene initial concentrations can be selected from the Type-2 and Type-3 PSO kinetic models in Table 5. In addition, the results showed that the equilibrium benzene degradation amount ( q e ) in PSO kinetic models were positively correlated with benzene concentration, and the rate constant K showed a trend of gradual increase with the pollutant concentration. The trend of the K value was opposite to the phenomenon observed by Zulfiqar et al. (Zulfiqar et al., 2019). This may be attributed to the increase in benzene concentration, which leads to higher collision probability between free radical particles and benzene molecules in the system, resulting in improved reaction rate (Shuchi et al., 2021).

In the IPD kinetic model, as the initial benzene concentration increased from 0.0625 mM, the value of diffusivity constant K id increased from 0.6002 to 9.334 mg/g min 0.5 (Table 5). This may be due to the stronger reaction driving force induced by the higher benzene concentration and thus affected the mass transfer rate in the reaction process (Cheung et al., 2007). However, it has been reported that the linear curves of negative boundary layer thickness ( C i ) and the linear curve of q t versus t 0.5 would not exceed zero at all concentrations (Abdelwahab and Amin, 2013). The IPD model was considered to be unable to reasonably explain the mechanism of the PS/nZVI oxidation technology due to abnormality in the fitted C i values (Figure 5D). Therefore, type II and type III PSO maybe the reasonable models to explain the benzene degradation process by the PS/nZVI oxidation system.