Simulation and Optimization of the CWPO Process by Combination of Aspen Plus and 6-Factor Doehlert Matrix: Towards Autothermal Operation

: This work aims to present an industrial perspective on Catalytic Wet Peroxide Oxidation (CWPO) technology. Herein, process simulation and experimental design have been coupled to study the optimal process conditions to ensure high-performance oxidation, minimum H 2 O 2 consumption and maximum energetic e ﬃ ciency in an industrial scale CWPO unit. The CWPO of phenol in the presence of carbon black catalysts was studied as a model process in the Aspen Plus ® v11 simulator. The kinetic model implemented, based on 30 kinetic equations with 11 organic compounds and H 2 O 2 involvement, was valid to describe the complex reaction network and to reproduce the experimental results. The computer experiments were designed on a six-factor Doehlert Matrix in order to describe the inﬂuence of the operating conditions (i.e., the di ﬀ erent process temperatures, inlet chemical oxygen demands, doses of H 2 O 2 and space time) on each selected output response (conversion, e ﬃ ciency of H 2 O 2 consumption and energetic e ﬃ ciency) by a quadratic model. The optimization of the WPO performance by a multi-criteria function highlighted the inlet chemical oxygen demand as the most inﬂuential operating condition. It needed to have values between 9.5 and 24 g L − 1 for autothermal operation to be sustained under mild operating conditions (reaction temperature: 93–130 ◦ C and pressure: 1–4 atm) and with a stoichiometric dose of H 2 O 2 .


Introduction
Advanced oxidation processes (AOPs) are technologies that involve the generation of hydroxyl radicals in sufficient quantities to mineralize recalcitrant pollutants to carbon dioxide, water and inorganics or, at least, to transform them into harmless products. Wastewater with relative low chemical oxygen demand (COD < 5 g L −1 ) is suitably treated under near ambient conditions by different types of oxidant i.e., H 2 O 2 , O 2 and O 3 . In addition, these processes can be photo and electro assisted. The treatment of highly loaded wastewater (COD > 20 g L −1 ) requires the employment of thermal technologies, namely wet air oxidation (WAO) or catalytic wet air oxidation (CWAO), that operate under high temperatures (T) and pressures (P), 130-300 • C and 5-200 atm, respectively [1]. The selection of each technology depends on many factors, such as the technical effectiveness using real streams, the manageability and applicability to the specific operations and the cost, among others [2][3][4].

Process Simulation Flowsheet
The simulation of the CWPO process was carried out with Aspen Plus ® v11. A flow and temperature diagram of the base case configuration is depicted in Figure 1. As can be seen, the wastewater, which contains the target pollutant (phenol in this case), enters the process (S-01) and mixes with a flow of H2O2 (S-02). The resulting stream (S-03) at TS-03 and a dose of H2O2 predefined with respect to the stoichiometric one for complete phenol mineralization are led to a heater (E-100) to increase the temperature to reaction conditions (TR). Then, the heated stream (S-04) is fed to the fixed-bed reactor (R-100), where the oxidation of phenol and its oxidized intermediates take place under isothermal conditions in the presence of a catalyst (carbon black, in this case). In the current process model, the reaction took place only in R-100. The CWPO reactor was modeled as an isothermal plug-flow reactor. The temperature of the liquid outlet stream of the reactor (S-05) was decreased in a cooler (E-101) to the temperature predefined for the process outlet stream (T-06). TS-06 was always 5 °C above TS-03 to maintain a minimum pinch point for the viability of the process heat exchanger network. The heat duty of the different unit operations was defined as H-01, H-02 and H-03, for the heater (E-100), cooler (E-101) and fixed-bed reactor (R-100), respectively ( Figure 1). The mass flows of the streams and dimensions of the reactor and the heaters were arbitrary, since only non-dimensional variables are taken into account in this work (see Section "Variables definition"). In addition, the pressure drop in the pipes and equipment was considered negligible, and the operating pressure was fixed at 4 bar, the same as that used to obtain the kinetic equations [40].

Reaction Pathway and Kinetic Model
The kinetic model used to simulate the CWPO of phenol upon a carbon black catalyst was experimentally obtained in a batch-reactor in the absence of mass transfer limitation and according to the reaction pathway previously proposed [40]. Briefly, phenol is oxidized to aromatic intermediates-namely, resorcinol (RS), catechol (CTL), hydroquinone (HQ) and p-benzoquinone (BQ)-in parallel reactions, and these cyclic species are broken down into different carboxylic acids, namely, maleic (MAL), fumaric (FUM), malonic (MLO), oxalic (OXA), acetic (ACE) and formic (FORM) acid. This reaction pathway comprises a set of 30   The mass flows of the streams and dimensions of the reactor and the heaters were arbitrary, since only non-dimensional variables are taken into account in this work (see Section "Variables definition"). In addition, the pressure drop in the pipes and equipment was considered negligible, and the operating pressure was fixed at 4 bar, the same as that used to obtain the kinetic equations [40].

Reaction Pathway and Kinetic Model
The kinetic model used to simulate the CWPO of phenol upon a carbon black catalyst was experimentally obtained in a batch-reactor in the absence of mass transfer limitation and according to the reaction pathway previously proposed [40]. Briefly, phenol is oxidized to aromatic intermediates-namely, resorcinol (RS), catechol (CTL), hydroquinone (HQ) and p-benzoquinone (BQ)-in parallel reactions, and these cyclic species are broken down into different carboxylic acids, namely, maleic (MAL), fumaric (FUM), malonic (MLO), oxalic (OXA), acetic (ACE) and formic (FORM) where C cat , C i and C H 2 O 2 refer to the concentrations of the catalyst, organic compound i and H 2 O 2 , respectively, while k i is the rate constant. The kinetic rate constant values were estimated from the fitted experimental data in the previous work of Diaz de Tuesta et al. [40]. The H 2 O 2 consumption rate was calculated from the production rate of each organic compound and the stoichiometric coefficient of the corresponding Equation (Table S1 of the Supplementary Material).

Components and Property Model Selection
The components considered in the simulation were H 2 O 2 , CO 2 , N 2 , O 2 , H 2 O, phenol and the oxidized species mentioned above. The UNIQUAC property model was selected to estimate the fluid properties in the process simulations. Most of the parameters were taken from the Aspen Property and NIST databases available in Aspen Plus. The unknown (water/organic compounds, CO 2 /O 2 and O 2 /N 2 ) binary interaction parameters were either estimated by the group contribution methods or fitted by the regression of experimental data using, respectively, the estimation and regression algorithms implemented in Aspen Plus. The experimental data of the corresponding vapor/gas-liquid equilibria were used to fit the unknown UNIQUAC binary parameters. They were retrieved from the literature [45][46][47][48][49] and the NIST database. The NIST database was accessed through the link implemented in Aspen Plus and supported by the NIST Thermo Data Engine. The binary interaction parameters used in this work are shown in Table S2 of the Supplementary Material. The gases considered in the simulations (CO 2 , O 2 and N 2 ) were defined as Henry components. The corresponding binary gas/water parameters for the Henry equation were obtained from the Aspen Property database. The heats of formation of all the components involved in the chemical reactions were taken from the Aspen Property database.

Variable Definition
The simulation of the CWPO process was run to explore the influence of six operating variables on six selected responses. The operating variables of the process studied were: T R , the reaction temperature ( • C). T input , the heater E-100 inlet temperature ( • C), equal to T S-01 = T S-02 = T S-03 as shown in Figure 1. T output , the cooler E-101 outlet temperature ( • C), equal to T S-06 as shown in Figure 1. ThOD feed , the COD of the wastewater (g·L −1 ). In this work, it is calculated from the concentration of phenol fed to the reactor then named as theoretical oxygen demand fed to the reactor. HP excess , the excess of H 2 O 2 with respect to the stoichiometric dose. τ, the space time in g cat mol phenol −1 h −1 , calculated as in Equation (2) by the ratio of the mass of the catalyst (W) to the inlet phenol molar flow rate (F).
The selected response variables were as follows: X j is the conversion of phenol, H 2 O 2 , theoretical oxygen demand (ThOD) or total organic carbon (TOC). It is calculated according to Equation (3), in terms of molar flow (for phenol and H 2 O 2 ) or mass flow (for ThOD and TOC): Catalysts 2020, 10, 548

of 17
The calculation of ThOD and TOC, g L −1 , for the inlet and outlet streams of the reactor is provided in η E is the energetic efficiency. It is calculated as the ratio of the process of "heat release" (through the cooler H-02 and reactor H-03) to "heat needed" (through the heater H-01). For that, a heat loss of 3% has been considered for heat exchangers and 10% for the chemical reactor [51]: Energetic efficiency values above 100% imply that the heat released by the CWPO reaction can be recovered to preheat the incoming feed stream. This means that the CWPO process is energetically self-sufficient or can be run under autothermal operation.

Design of Experiments
The study of the influence of the selected operating variables T R , T input , T output , ThOD feed , HP excess and τ on the responses X phenol , X ThOD , X H 2 O 2 , X TOC , η H 2 O 2 and η E in the CWPO process, with the minimum number of simulation runs, was carried out based on a six-factor Doehlert Matrix.
The resulting numerical experiments, or cases, to simulate by Aspen Plus numbered 43. Table 1 shows these experiments with the real and code values for each operating variable. The different ranges of interest to assign the real values were the following: T R from 50 to 130 • C, taking into account that the oxidation rate is almost negligible below 50 • C [40]; ThOD feed from 2.38 to 23.8 g·L −1 , corresponding to an inlet phenol concentration equal to 1−10 g L −1 , an appropriate pollutant load for the CWPO processes [8,36]; T output from 30 to 60 • C, taking into account the values of the pitch point above 5 • C-the minimum value is 8.7 • C in Case 35 and the maximum, 41.2 • C in Case 41. Additionally, the Doehlert Matrix allows the simulation of this range of values without temperature crossing, thus avoiding inconsistent behavior scenarios such as T output values above T R . The HP excess varied from 50% to 150% respect to the stoichiometric dose. The τ was between 94 and 235 g cat ·h·mol phenol −1 , in the range of values reported for flow reactors [15][16][17][18][19][20][21][22]. T input ranged from 15 to 25 • C, which reflects the worst scenario. The function selected to describe the dependence between each response and the operating variables was a second-order polynomial equation: where y represents the predicted response; x i and x j , the codes of the operating variables (provided in Table 1 as x 1 = T R , x 2 = ThOD feed , x 3 = T output , x 4 = HP excess , x 5 = τ and x 6 = T input ); and β, the regression coefficients where β 0 is the intercept, β i is the coefficient for linear effects, β ii is that for quadratic effects and β ij is that for two-factor interaction effects. The Solver tool in Microsoft ® Excel, which employs the Generalized Reduced Gradient algorithm, was used to estimate the coefficients by multiple non-linear regression. The analysis of variance (ANOVA) was used to evaluate both the coefficients' statistical significance and the fitting of the function. The objective functions were maximized or minimized, as appropriate, using as restrictions the minimum and maximum values of the operating variables provided in Table 1. The optimization was run under several interactions and with more than 100 starting points (arbitrarily chosen).

Kinetic Model Validation
Prior to the simulation of the CWPO process, it is necessary to confirm that the evolution of the reactants (i.e., phenol and H 2 O 2 ) and the oxidized intermediates can be successfully computed, and that the results obtained reproduced the experimental ones. Based on the elementary reaction set introduced in the simulator (Table S1), the empirical kinetic equation was defined for each organic species, shown in Equation (1), with the kinetic constants estimated for the CWPO with a carbon black catalyst [40] and the thermodynamic data selected (Table S2 of the Supplementary Material); calculations were performed in a batch reactor with Aspen Plus under the following selected operating conditions: C Phenol,0 = 1 g·L −1 , C H 2 O 2 ,0 = 5 g·L −1 , pH 0 = 6, W = 4 g cat and T = 90 • C and at different reaction times in order to mimic the experimental conditions.
The results obtained in the simulations, in terms of the temporal evolution of the molar fraction of phenol and the aromatic intermediates (BQ, CTL, HQ and RS), are provided in Figure 2a. The profiles for the organic acids, also calculated, are not shown for simplicity. Besides, the H 2 O 2 molar fraction calculated from the production rate of each organic compound and the stoichiometric coefficient of the corresponding equation (Table S1) is provided in Figure 2b. As can be seen in Figure 2c, the molar fractions calculated by the simulation agree with the experimental data obtained in our previous work [40]. The kinetic model implemented in the simulator is valid to describe the complex reaction network of the phenol CWPO in the presence of carbon black catalysts. Note that the kinetic model is particular for each catalyst, since the specificity of each active site can lead to a different reaction networks and rate equations.
that the results obtained reproduced the experimental ones. Based on the elementary reaction set introduced in the simulator (Table S1), the empirical kinetic equation was defined for each organic species, shown in Equation (1), with the kinetic constants estimated for the CWPO with a carbon black catalyst [40] and the thermodynamic data selected (Table S2 of the Supplementary Material); calculations were performed in a batch reactor with Aspen Plus under the following selected operating conditions: CPhenol,0 = 1 g·L −1 , CH 2 O 2 ,0 = 5 g·L −1 , pH0 = 6, W = 4 gcat and T = 90 °C and at different reaction times in order to mimic the experimental conditions.
The results obtained in the simulations, in terms of the temporal evolution of the molar fraction of phenol and the aromatic intermediates (BQ, CTL, HQ and RS), are provided in Figure 2a. The profiles for the organic acids, also calculated, are not shown for simplicity. Besides, the H2O2 molar fraction calculated from the production rate of each organic compound and the stoichiometric coefficient of the corresponding equation (Table S1) is provided in Figure 2b. As can be seen, the molar fractions calculated by the simulation agree with the experimental data obtained in our previous work [40]. The kinetic model implemented in the simulator is valid to describe the complex reaction network of the phenol CWPO in the presence of carbon black catalysts. Note that the kinetic model is particular for each catalyst, since the specificity of each active site can lead to a different reaction networks and rate equations.

CWPO Performance
The values of the response variables (i.e., Xphenol, XThOD, XTOC, XH2O2, ηH2O2 and ηE) obtained in the 43 simulation cases conducted under the operating conditions defined by the Doehlert Matrix (Table

CWPO Performance
The values of the response variables (i.e., X phenol , X ThOD , X TOC , X H 2 O 2 , η H 2 O 2 and η E ) obtained in the 43 simulation cases conducted under the operating conditions defined by the Doehlert Matrix (Table 1) are given in Table 2. As observed, total pollutant removal (X phenol > 99%) is only achieved in Cases 2, 4, 14, 31 and 32, in which X ThOD and X TOC are also high, above 75% and 71%, respectively, as is η H 2 O 2 , between 69% and 79%. There is always remaining ThOD and TOC in the reactor outlet. In fact, the maximum conversions are obtained at the highest reaction temperature of 110 • C, in Cases 4 and 14, where both X ThOD~7 8% and X TOC~7 4% are achieved. Noteworthily, the self-sustaining thermal condition (η E > 100%) is achievable under a wide range of operating conditions (in 24 scenarios of the 43 simulated, see Table 2), including the whole ranges of values explored for T output (30-60 • C), T input (15-25 • C), HP excess (50%-100%) and τ (94 and 235 g cat ·h·mol phenol −1 ). This is an important aspect because in these 24 scenarios, the CWPO unit can be operated without any external heat supply under steady state operation. The ThOD feed seems to be the most influential variable affecting the energetic efficiency (see Figure S1 of the Supplementary Material). Additionally, it is possible to observe a linear relationship between η E and ThOD feed , regardless of the values of the rest of the operating conditions-T R , T output , HP excess , τ and T input (see Figure S1d of the Supplementary Material). This can be explained because the more the organic matter that is converted in the reactor (by the oxidation of phenol and its intermediates), the more heat is generated (the same tendency was observed with respect to X TOC ). However, ThOD feed should not be the only variable taken into account, as exemplified by Cases 36 and 30. In the former, 5.4 g L −1 of ThOD was removed in the reactor (ThOD feed = 13 g/L and X ThOD = 42%) and the energy efficiency was 101%, while in the latter, 8.4 g L −1 of ThOD was removed in the reactor (ThOD feed = 13 g·L −1 and X ThOD = 64.7%) while η E was lower than 100%. Elucidating the interactive effects of the operating variables on the responses will necessitate assessing the mathematical adjustment of the data to the quadratic model in Equation (6), as addressed below. Finally, the minimum X phenol , X ThOD and X TOC values that guarantee the self-sustaining thermal CWPO condition (η E > 100%) are 75, 42 and 37%, respectively (Case 36 of Table 2). In this case, an energy efficiency of 101% is achieved at ThOD feed = 13 g/L, HP excess = 60%, τ = 153 g cat ·h·mol phenol −1 ,

Mathematical Modeling and Response Surface Plots
The ANOVA results of quadratic model for X phenol , X ThOD , X TOC , X H 2 O 2 , η H 2 O 2 and η E are summarized in Table 3, along with the estimated coefficient values. The F-distribution and p-values for each estimated coefficient are collected in order to compute their statistical significance (coefficient excluded if p-value > 0.30). As observed, the F of Fisher is higher than the p-value for all the regressions, so the desired model significance is fulfilled. In fact, the parity plot depicted in Figure 3 shows good agreement.  a higher concentration of radical species in the reactor medium increases the organic compound oxidation rates but also favors the conversion of the radical scavenging recombination into H2O2. The most influential variable for the conversions is ThODfeed (β2 > βi≠2), which is in agreement with the previous discussion.
According to this, a high reaction temperature improves the efficient consumption of H2O2 in CWPO, as has been already reported for the Fenton process [52]. To a lesser extent, a high concentration of ThODfeed also favors the efficient use of the oxidant. Note that doses of H2O2 above the stoichiometric ones, HPexcess > 100%, do not induce a significant detriment to the efficiency. In the case of ηE, the significant operating conditions are Tinput, Toutput and ThODfeed-in particular, ThODfeed (β2 = 63.2 vs. |βi≠2| = 8-25.4). Expectedly, a higher initial pollutant concentration in the feed stream leads to higher conversions, and thus, more heat is released in the reactor. It is also noteworthy  In general, the conversion response equations are similar since all the conversions (i.e., X phenol , X ThOD , X TOC and X H 2 O 2 ) take high values upon increasing the T R (β 1 = 15.2-22.1), ThOD feed (β 2 = 35.1-40.8) and τ (−β 5 = 10.2-12.5). Interestingly, X phenol , X ThOD and X TOC also increase upon increasing HP excess (β 4 = 16.6-18.4), whereas the X H 2 O 2 is decreased by increasing HP excess (β 4 = −13.1), which indicates that a higher concentration of radical species in the reactor medium increases the organic compound oxidation rates but also favors the conversion of the radical scavenging recombination into H 2 O 2 . The most influential variable for the conversions is ThOD feed (β 2 > β i 2 ), which is in agreement with the previous discussion.
Regarding η H 2 O 2 , T R is the variable with the greatest influence (β 1 > β i 1 ) followed by ThOD feed . According to this, a high reaction temperature improves the efficient consumption of H 2 O 2 in CWPO, as has been already reported for the Fenton process [52]. To a lesser extent, a high concentration of ThOD feed also favors the efficient use of the oxidant. Note that doses of H 2 O 2 above the stoichiometric ones, HP excess > 100%, do not induce a significant detriment to the efficiency.
In the case of η E , the significant operating conditions are T input , T output and ThOD feed -in particular, ThOD feed (β 2 = 63.2 vs. |β i 2 | = 8-25.4). Expectedly, a higher initial pollutant concentration in the feed stream leads to higher conversions, and thus, more heat is released in the reactor. It is also noteworthy that while T input has a positive effect (β 6 = 10.0), T output has a negative one (β 3 = −25.4); this means that η E is increased by decreasing the heat recovery from the output stream. For this response, the effect of this temperature (T output ) is found to be more significant than the effect of the T R , ThOD feed , HP excess or τ (β 6 > β 1 , β 3 , β 4 , β 5 ).
All the fitted response variables, i.e., X phenol , X ThOD , X TOC , X H 2 O 2 , η H 2 O 2 and η E , show significant second-order (β ii ) and interaction (β ij ) coefficients, so it will be important to consider all the operating conditions to seek optimal performance. To gain an insight into this, the response surface plots and 2-D contour plots showing the effects of T R and ThOD feed and also, HP excess and ThOD feed on the responses X phenol , X ThOD , η H 2 O 2 and η E are shown in Figures 4 and 5, respectively.
As can be seen, for X phenol (A) and X ThOD (B), the increase in T R and HP excess up to 100 • C and 100% of the stoichiometric dose, respectively, enhances the conversions, whereas ThOD feed values above 16 g·L −1 have a clearly negative impact. This means that X phenol and X ThOD can be optimized to achieve their maximum values (100 and 78%, respectively) inside the operation window considered.
Regarding the η H 2 O 2 , its variation is higher with T R than ThOD feed and HP excess , as expected, but the three variables must be considered for the optimization of η H 2 O 2 . It is remarkable that the thermal self-sustaining conditions (η E ≥ 100%) can be achieved across a wide range of η H 2 O 2 values (47%-81%) (see Figures 4c and 5c). Lastly, η E ≥ 100% are obtained at T R from 80 to 100 • C and ThOD feed values above 9 g·L −1 (Figure 4d). The effect of the HP excess is not as significant as that of ThOD feed (Figure 5d) for this response.
Catalysts 2020, 10, x FOR PEER REVIEW 11 of 17 that while Tinput has a positive effect (β6 = 10.0), Toutput has a negative one (β3 = −25.4); this means that ηE is increased by decreasing the heat recovery from the output stream. For this response, the effect of this temperature (Toutput) is found to be more significant than the effect of the TR, ThODfeed, HPexcess or τ (β6 > β1, β3, β4, β5). All the fitted response variables, i.e., Xphenol, XThOD, XTOC, XH 2 O 2 , ηH 2 O 2 and ηE, show significant second-order (βii) and interaction (βij) coefficients, so it will be important to consider all the operating conditions to seek optimal performance. To gain an insight into this, the response surface plots and 2-D contour plots showing the effects of TR and ThODfeed and also, HPexcess and ThODfeed on the responses Xphenol, XThOD, ηH 2 O 2 and ηE are shown in Figures 4 and 5, respectively. (%) As can be seen, for Xphenol (A) and XThOD (B), the increase in TR and HPexcess up to 100 °C and 100% of the stoichiometric dose, respectively, enhances the conversions, whereas ThODfeed values above 16 g·L -1 have a clearly negative impact. This means that Xphenol and XThOD can be optimized to achieve their maximum values (100 and 78%, respectively) inside the operation window considered. Regarding the ηH 2 O 2 , its variation is higher with TR than ThODfeed and HPexcess, as expected, but the three variables must be considered for the optimization of ηH 2 O 2 . It is remarkable that the thermal selfsustaining conditions (ηE ≥ 100%) can be achieved across a wide range of ηH 2 O 2 values (47%-81%) (see Figure 4c and 5c). Lastly, ηE ≥ 100% are obtained at TR from 80 to 100 °C and ThODfeed values above 9 g·L −1 (Figure 4d). The effect of the HPexcess is not as significant as that of ThODfeed (Figure 5d) for this response.

Multi-Criteria Optimization of CWPO Performance
The operating conditions necessary to obtain the highest Xphenol, XThOD, XTOC, XH 2 O 2 , ηH 2 O 2 or ηE upon CWPO in the presence of carbon black catalysts, derived from the optimization of the quadratic model selected in Equation (6), for each response are provided in Table 4.

Multi-Criteria Optimization of CWPO Performance
The operating conditions necessary to obtain the highest X phenol , X ThOD , X TOC , X H 2 O 2 , η H 2 O 2 or η E upon CWPO in the presence of carbon black catalysts, derived from the optimization of the quadratic model selected in Equation (6), for each response are provided in Table 4. As can be seen, the space time, τ, required to achieve this best individual performance for each response needs to be high, usually 235 g cat ·h·mol phenol −1 . In order to assure the maximum high degrees of oxidation in terms of phenol, oxygen demand and TOC, only wastewater with high ThOD feed (>16.2 g·L −1 ) could be suitably treated under high reaction temperatures (100-130 • C) and an excess of H 2 O 2 dose. However, this scenario is far from that in which complete and efficient H 2 O 2 consumption occurs and, consequently, far from being a cost-effective treatment [2,36]. The wastewater with the maximum chemical oxygen demand selected (ThOD feed = 23.8 g·L −1 ) and its treatment under an excess of H 2 O 2 guarantee the production of energy upon CWPO (η E > 200%) while the reaction temperature can be set to 65 • C, and T input and T output to 25 and 30 • C, respectively (Table 4).
Considering the practical application of the treatment, the CWPO performance has been optimized in order to achieve the maximum oxidation degree, in terms of X ThOD , with the most efficient consumption of H 2 O 2 and most favorable energetic conditions. For that, a multi-objective optimization has been carried out using linear scalarization to find the optimal values for the operating conditions (T R , HP excess , τ and T output ) for a given wastewater (T input and ThOD feed ) that assure the greatest values of the selected objectives (X ThOD , η H 2 O 2 and η E ): where y opt is the function that allows increasing the three responses X ThOD , η H 2 O 2 and η E with the same weight in Equation (7). The optimal operating conditions for T R , T output , HP excess and τ in the 10 cases of multi-criteria optimization calculated by setting T input = 20 • C and ThOD feed between 2.38-23.8 g·L −1 are summarized in Table 5. Besides, the values of the six response variables-i.e., X phenol , X ThOD , X TOC , X H 2 O 2 , η H 2 O 2 and η E -estimated by simulation with Aspen Plus and by the response quadratic models (in brackets) are included.
According to the above results, the CWPO process should be considered as a treatment able to abate recalcitrant compounds such as phenol and to increase the detoxification of the wastewater (due to the oxidation of recalcitrant compounds into organic acids) before its disposal via a biological process. The CWPO process could be successful in treating a wide range of industrial wastewaters as a one-step treatment and then oriented to increase the quality of the treated effluent for water reuse or safe discharge.

Conclusions
It has been demonstrated that CWPO aims to fill the gap between the AOPs under near ambient conditions and CWAO processes. The advantages of CWPO include the treatment of mildly loaded wastewater, with a COD demand between 9.5 and 24 g L −1 , under mild operating conditions (T R = 93-130 • C, P = 1−4 atm and a stoichiometric dose of H 2 O 2 ) in an autothermal fixed bed reactor. In addition, there is the possibility of energy recovery, to a greater or lesser extent depending on the inlet COD and the reaction temperature.
This study demonstrates that: • The COD removal achieved in CWPO units is increased by increasing the reaction temperature and dose of hydrogen peroxide. Nevertheless, the best conditions lead to 78% COD conversion due to the refractoriness of the organic acid species.

•
Working at high temperature favors the efficient consumption of H 2 O 2 in the oxidation process.

•
Autothermal operation is possible if a high load of COD (or TOC) is employed and high conversions are achieved, which requires setting up the process with the optimal operating conditions.
Clearly, the CWPO process is a promising technology for reducing the impact of non-biodegradable hazardous substances on the environment.  Table S2: Binary coefficients (aij, aji, bij, bji) estimated for the implementation of the UNIQUAC thermodynamic model in the simulator software, Table S3: Theoretical oxygen demand (ThOD) and total organic carbon (TOC) of each compound (i).