RSM experimental design and hydrocarbon degradation in planted mesocosms
First, the test values of the five variables nutrients (A), surfactant (B), aeration (C), hydrocarbon content (D), and hydraulic retention time (E) were chosen at three levels [low (− 1), central (0), and high (+ 1)] based on previous studies (Table 1)11,20,21,22,23. Then, CCD was used to generate the experimental design matrix. CCD was favored over a Box Behnken design (BBD) because it offers more axial design points compared to the BBD while being suitable for testing five variables34. Furthermore, CCD is better at extreme conditions and gives better results for quadratic models35. In this study, the matrix consisted of 32 factorial points, 10 axial points, and 8 central points, resulting in a total of 50 experimental setups (Table 2).
The 50 different setups were established in triplicates as 3-L mesocosms with hydroponically grown common reed (Phragmites australis) (Fig. 1). Table 2 shows the results of hydrocarbon removal (% concentration reduction), COD reduction (%), and growth of plant biomass (g) in the setups. The highest hydrocarbon removal (89%) occurred with A: 14 mg L−1 nitrogen and 1.9 mg L−1 phosphorus; B: 0.005% (w/v) of sodium dodecyl sulphate as surfactant; C: 1 L of air min−1; D: 0.75% hydrocarbon content; and E: 24 days (setup #46). The lowest hydrocarbon removal among all 50 setups was 6% (setup # 38 = HC: 0.5 mg/L, surfactant: 0%, aeration: 0 L/min, nutrients ratio: 0, and retention time: 8 days), and the lowest removal with a hydraulic retention time of 24 days was 19% (setup # 18 = HC: 0.5 mg/L, surfactant: 0%, aeration: 2 L/min, and nutrients ratio: 2), which was a substantial difference among these two setups.


Preparation of mesocosms with P. australis for the optimization of FTW operational parameters.
Development and validation of response surface models
Next, RSM was applied for the mathematical model building of the experimental data obtained with the hydroponic systems. Multiple regression analysis was used to test for linear (A, B, C, D, E), quadratic (A2, B2, C2, D2, E2) and interactive (AB, AC, AD, AE, BC, BD, BE, CD, CE, DE) effects of all variables. The following polynomial quadratic equations fitted best to the experimental response data for hydrocarbon degradation, COD reduction, and increase in plant biomass (Eqs. 1–3).
$${text{Y}}_{{{text{hydrocarbon reduction (% )}}}} = , + { 8}.{35 }{-} , 0.{2}0{text{5 A }} + , 0.{text{689 B }} + , 0.{5}0{text{6 C }} + , 0.{text{417 D }} + { 1}.{text{53 E }}{-} , 0.{text{186 BC}}{-}0.{25}0{text{ BD }}{-} , 0.{text{229 CD }}{-}{ 1}.{text{13 C}}^{{2}} {-}{ 1}.{text{12 D}}^{{2}} .$$
(1)
$${text{Y}}_{{{text{COD reduction }}(% )}} = , + { 78}.{62 }{-}{ 2}.{text{79 A }} + { 8}.{text{32 B }} + { 6}.{text{15 C }} + { 4}.{text{56 D }} + {19}.{text{29 E }}{-}{ 2}.{text{78 BD }} + { 2}.{text{59 DE }}{-}{ 13}.{text{49 C}}^{{2}} {-}{ 15}.{text{49 D}}^{{2}} .$$
(2)
$${text{Y}}_{{{text{plant biomass }}({text{g}})}} = , + { 34}.{75 }{-}{ 1}.{text{12 A }} + { 2}.{text{12 B }} + , 0.{text{9118 C }} + { 1}.{text{26 D }} + { 3}.{text{41 E }} + , 0.{text{937 BE }} + , 0.{text{875 CE }} + , 0.{text{875 DE }}{-}{ 3}.{text{49 A}}^{{2}} {-}{ 6}.{text{99 C}}^{{2}} {-}{ 6}.{text{99 D}}^{{2}} + { 8}.{text{51 E}}^{{2}} .$$
(3)
where Y is the response value, A stands for nutrient concentration, B for surfactant concentration, C for aeration, D for hydrocarbon content, and E for retention time; AB, AC, AD, AE, BC, BD, BE, CD, CE, and DE are the interaction effects; A2, B2, C2, D2, and E2 represent square effects. The negative (−) and positive (+) signs of regression coefficients showed that there were antagonistic and synergistic effects of the variables. Insignificant terms with p > 0.05 were removed from the three models.
All five variables possessed the same linear significant terms A, B, C, D, and E and quadratic terms C2 and D2. The interaction terms BC, BD, and CD were significant for hydrocarbon reduction, BD and DE were significant for COD reduction and BE was the only significant interaction term for the production of plant biomass.
An analysis of variance (ANOVA) confirmed the adequacy of the quadratic models for the three responses with p-values < 0.0001 (Table 3). Precisely, hydrocarbons attenuation, COD reduction, and gain in plant biomass were tested by fitting quadratic models in RSM. This approach describes the mathematical relationship between each term in the model and response. The coefficient of determination (R2) was 0.95 for the attenuation of the hydrocarbon concentration (Fig. 2a). For COD reduction and for increase in plant biomass it was R2 = 0.93 and R2 = 0.88, respectively (Fig. 2b,c). The independent variables accounted for 96% of the variability. Furthermore, there were strong relations of surfactant, aeration, and nutrients with R2 values of 0.95, 0.939, and 0.883, respectively (Table 2). The goodness-of-fit of the regression equation was confirmed by the high value of the adjusted determination coefficient (R2adj = 0.938). This high value showed that the selected factors and their values constitute a very good representation of the main processes that influence the hydrocarbon treatment efficiency of the FTW systems.


(a–c) Correlation between the actual and predicted (a) hydrocarbon reduction, (b) COD reduction, and (c) increase in plant biomass.
Model analysis via 2D contour graphs and 3-D surface plots
To visualize the relationships between the experimental variables and the corresponding responses, we used RSM to draw three-dimensional response surface graphs and contour plots as their two-dimensional projections (Fig. 3). In this approach, the significance of mutual effects of the experimental variables is represented by the curvature of the response surface and contour lines. Saddle and ridge-shaped 3-D graphs [(inverse) hyperbolic contour plots] exhibit a significant mutual interaction of experimental variables, while a dome-shaped 3-D graph (circular contour plot) represents a non-significant interaction. Here, the effect of the experimental factors was investigated by varying two factors over the experimental range while keeping the other three variables constant.


Response surface graphs for hydrocarbon reduction, COD reduction, and increase in plant biomass.
Figure 3 illustrates the effect of the variables on hydrocarbon decrease. The surfactant concentration and aeration significantly affected hydrocarbon reduction by varying levels of both variables. The 3-D diagram displays that hydrocarbon attenuation increases with increasing surfactant concentration whereas an increase in the level of aeration helps to decrease hydrocarbon concentration up to an optimum point (Fig. 3a). Similar results were found for the interaction between surfactant and nutrients (Fig. 3b). The ridge shape of the 3-D graph is showing a significant interaction of the variables. Higher surfactant concentrations produced a positive effect on the degradation of hydrocarbons while higher nutrient concentrations resulted in an increase in hydrocarbon attenuation to an optimum point, after which further nutrient increase caused a negative effect on the response in the model. In Fig. 3c the dome surface of the response plot shows that the interaction between nutrients and aeration is non-significant. The 3-D surface plot indicates that both high and low levels of nutrients and aeration did not have a statistically significant effect on hydrocarbon degradation.
The interactive effect of the experimental variable on COD reduction was also determined using 3D plots of RSM. The 3-D graph shows that the interaction between the variables is significant. Higher levels of surfactant had a positive effect to decrease COD in the water whereas after an optimum level a further increase in nutrients has a negative effect on the process (Fig. 3d). As expected, COD was most effectively reduced at the highest level of retention time, as shown in Fig. 3d. The effect of the variables on the growth of plant biomass was also demonstrated by the design expert. The 3-D graphs in Fig. 3e,f show that an increase in retention time increases the plant biomass significantly, while the various levels of surfactant and aeration have static or limited effects on plant biomass. A similar trend was observed for retention time, nutrients, and aeration (data not shown).
Optimization of experimental conditions for hydrocarbon degradation
Then, RSM was used to predict the optimal values of the variables namely nutrients, surfactant, aeration, hydrocarbon content with a hydraulic retention time of 24 days to maximal attenuation of the hydrocarbon concentration. The optimized values of variables predicted by the desirability function method of RSM were found to be a hydrocarbon content of 0.758%, a surfactant concentration of 0.006%, aeration of 1.178 L of air min−1, and a nutrient ratio of 1.20, resulting in a predicted value of hydrocarbon degradation of 98% (Fig. 4). Then we carried out another experimental test at the 3-L scale with the optimized values predicted by RSM. Attenuation of the hydrocarbon concentration of 95% was achieved in the FTW setup with the RSM-optimized operational parameters. Thus, the experimentally observed response values agree again very well with the theoretical values assumed by the model, showing the precision and accuracy of the RSM approach.


Desirability ramps for numerical optimization of hydrocarbon reduction, COD reduction and plant biomass.
The benefit of RSM for improved hydrocarbon degradation in FTW
In general, an RSM model can be used to predict what will happen under different conditions, but it cannot explain the mechanism of the process34. Nevertheless, the goodness of fit between the predicted and experimental values can indicate whether all important parameters have been accounted for in the model, and thus whether the underlying conceptual process framework is close to reality. As reported above, the adjusted determination coefficient for the model equations in this study were R2adj = 0.938 with probability values of p < 10–4, demonstrating the significance of the model to predict the responses and thus fostering a rational and cost-effective improvement of FTW-based water treatment of oil-contaminated water at field scale. It is important that the two operational parameters retention time and level of aeration could be successfully modeled with RSM, as these are prime parameters in process engineering. These parameters of the system can be more readily adjusted to achieve the desired treatment efficiency at given operational costs. It is also important to note that there was an optimal aeration level. To consider this finding may limit costs at field-scale applications.
There are two potential limitations of the present study for translating its results to full-scale systems. First, the present investigation was carried out in batch mode. Several studies with FTW at scales ranging from laboratory to field scale have shown that results gained at smaller scales are essentially valid for the field scale, however, it is not a given that this is always the case. Secondly, long-term effects were not investigated in this study. The removal of hydrocarbons during the continuous operation of FTWs will have to be investigated in future work. Aspects of the FTW such as vitality of plants, dimensions of root network, i.e., a volume ratio of root network to free water, will change over time and may affect treatment performance.
Cost-effectiveness ratio (CER)
In this study, CER was estimated yearly. At first, the total present value cost (pvc) for a single 1000 L wetland system was calculated as Eq. (4).
$${text{pvc }}left( {$ {581}} right) , = {text{ pvc}}_{{{text{ic}}}} left( {$ {456}} right) , + {text{ pvc}}_{{{text{om}}}} ; left( {$ {125}} right).$$
(4)
Then, CERtotal for 12 months of operation was calculated by dividing pvc by the volume of water receiving treatment, multiplying the number of required treatments (n) (Eq. 5).
$$CER, left(totalright)=frac{PV, of, total, costs left($581right)}{ volume ,of, water, receiving, treatment ; left(12,000, {text{L}}right)}times text{ n},$$
(5)
where n is a factor representing the number of times the system has been operated.
This indicated that, by following RSM, we can reduce the treatment cost up to $0.048n per liter of total wastewater receiving treatment.
RSM was successfully applied to optimize the abiotic variables nutrients concentrations, surfactant addition, aeration, and retention time for the attenuation of hydrocarbons from oil-contaminated water in a mesocosm-scale FTW experiment. The optimum values of the operational parameters were at a crude oil concentration of 0.758%, aeration 1.178 L of air min−1, a surfactant concentration of 0.006%; a nutrient ratio of 1.20; and retention time of 23.6 days for maximum hydrocarbons removal from the water, which resulted in a predicted and experimental attenuation of 98% and 95%, respectively. The performance of the system mainly depended on the retention time, but the initial oil concentration, surfactant concentration, nutrient ratio, and aeration rate also affected the removal of hydrocarbons from the water. Effect of salinity in crude oil wastewater treatment is nevertheless crucial, which may be included in the RSM design for future studies. Also, the results of RSM efficacy should be validated at pilot- and/or operational-scale for field-oriented conclusions. Thus, this study shows that the use of RSM is promising for reducing the costs of field-scale operation of FTW for hydrocarbon attenuation at oil processing sites.

