Operating Conditions Optimization via the Taguchi Method to Remove Colloidal Substances from Recycled Paper and Cardboard Production Wastewater

Optimization of the ultrafiltration (UF) process to remove colloidal substances from a paper mill’s treated effluent was investigated in this study. The effects of four operating parameters in a UF system (transmembrane pressure (TMP), cross-flow velocity (CFV), temperature and molecular weight cut-off (MWCO)) on the average permeate flux (Jv), organic matter chemical oxygen demand (COD) rejection rate and the cumulative flux decline (SFD), was investigated by robust experimental design using the Taguchi method. Analysis of variance (ANOVA) for an L9 orthogonal array were used to determine the significance of the individual factors, that is to say, to determine which factor has more and which less influence over the UF response variables. Analysis of the percentage contribution (P%) indicated that the TMP and MWCO have the greatest contribution to the average permeate flux and SFD. In the case of the COD rejection rate, the results showed that MWCO has the highest contribution followed by CFV. The Taguchi method and the utility concept were employed to optimize the multiple response variables. The optimal conditions were found to be 2.0 bar of transmembrane pressure, 1.041 m/s of the cross-flow velocity, 15 °C of the temperature, and 100 kDa MWCO. The validation experiments under the optimal conditions achieved Jv, COD rejection rate and SFD results of 81.15 L·m−2·h−1, 43.90% and 6.01, respectively. Additionally, SST and turbidity decreased by about 99% and 99.5%, respectively, and reduction in particle size from around 458–1281 nm to 12.71–24.36 nm was achieved. The field-emission scanning electron microscopy images under optimal conditions showed that membrane fouling takes place at the highest rate in the first 30 min of UF. The results demonstrate the validity of the approach of using the Taguchi method and utility concept to obtain the optimal membrane conditions for the wastewater treatment using a reduced number of experiments.


Introduction
The pulp and paper (P & P) industry is ranked as the world's third largest consumer of fresh water [1] and an important producer of wastewater with different organic and inorganic contaminants. Depending on the type of processes used in paper manufacture, the integration between production and environmental protection is one of the key topics in the paper industry.
According to the Confederation of European Paper Industries (CEPI) [2], Europe is the second largest producer of paper and paperboard with 22.7% (91.39 million tons) of world production, of multiple response optimization, because the process performance is often evaluated using several quality characteristics (responses). In this case, the Taguchi method and utility concept are useful tools for optimizing operating parameters in multiple characteristics responses [34].
The statistical analysis of variance (ANOVA) can be used to provide information on whether the operating parameters (factors) are statistically significant or not, as well as to identify the influence of individual factors and establish the relationships between the factors and operating conditions. In this analysis, the p-value index is used to know which operating parameters have a significant effect on the response variables. This information is complemented with the F-test which helps to identify whether or not each factor is significant at the selected confidence level [35,36]. In this study, ANOVA was also used to analyze the experimental results.
The aim of this work was to determinate the effect of operating conditions such as transmembrane pressure (TMP), cross-flow velocity (CFV), temperature and molecular weight cut-off (MWCO) on the average permeate flux, chemical oxygen demand (COD) rejection rate and cumulative flux decline (SFD), in addition to determining the optimal conditions for the given sets of values and to find the best response variables by using Taguchi experimental design and the utility concept. The results of this study may be used as a guideline when operating UF systems under the best conditions in a wastewater treatment plant (WWTP) in a papermaking factory. The filtration results and analysis of the experimental data presented and discussed in this study were carried out by using ANOVA to find the significance of the controlling factors and optimized using the Taguchi method to find the optimal operating conditions. A standard L 9 orthogonal array was selected for experimental planning with four factors and three levels for each factor.

Paperboard Mill Treated Effluent Feedstock
The paperboard mill treated effluent (PMTE) used in this work came from a secondary clarifier effluent from a WWTP in a papermaking factory located in the south of the Valencian autonomous region in Spain. In order to prevent early membrane fouling, remove large suspended solids, and reduce initial turbidity and COD in the PMTE, the raw feed solution was pre-filtered by conventional filtration (low-pressure pump at around 1 bar) with a Cintropur ®® NW 50 filter element, and centrifugal propeller and filter cloths with a 5 µm nominal pore size. The significant characteristics of the PMTE samples are listed in Table 1. Table 1. Average compositions of the ultrafiltration (UF) feed solution (prefiltered solution from a secondary clarifier effluent from a wastewater treatment plant (WWTP)).
The experiments were performed in a typical UF pilot plant with a flat-sheet membrane module (Rhône-Poulenc, Lyon, France), that allowed working with two membranes with similar or different MWCO (depending on the experiment being carried out). The effective area for each membrane in Membranes 2020, 10,170 4 of 22 the module was 154.8 cm 2 . The details of the experimental set-up have been described previously by Sousa et al. [21]. The pilot plant had a data acquisition system (temperature, module input and output pressure) from LabVIEW ® (National Instruments, Austin, TX, USA). The permeate was collected during the filtration in a beaker placed on an electronic balance connected to a computer in order to continuously register the weighting data. This data was then automatically logged every 30 s and subsequently used to calculate permeate flux through the membranes.

Analytical Methods
The PMTE used as the feed solution and the UF permeate samples were analyzed according to the methods described below. The suspended solids analyses were carried out in accordance with the standard methods [37]. Turbidity was measured using a Dinko 112 turbidimeter (ASTM D1889, Barcelona, Spain). Conductivity was measured using a WTW level 3 conductivity device (ASTM D1125-14, Weilheim in Oberbayern, Germany). COD and total nitrogen in the effluent were analyzed using a Merck photometer (Merck KGaA, Darmstadt, Germany) and a Merck TR-300 thermoreactor (Merck KGaA, Darmstadt, Germany) in accordance with the standard methods.

Field-Emission Scanning Electron Microscopy (FESEM)
Field-emission scanning electron microscopy (FESEM) measurements were used to provide information on the fouling that formed on the membranes. The surface and cross-section morphologies of the fresh and fouled membranes were observed by field-emission scanning electron microscope, (ZEISS ULTRA 55, Oxford Instruments, Berkshire, UK), operated with a voltage of 200 kV and an accelerating voltage of 0.02-5 kV. Before analysis, the dried membrane samples were attached to double-sided adhesive carbon tape on an aluminum holder, and subsequently coated with a thin layer of gold prior to analysis.

Membrane Characterization
Before the UF runs, permeability experiments were carried out to determine the intrinsic membrane resistance (R m ). Distilled water was used as the feed solution and measurements were taken for 1.0, 1.5, 2.0, and 3.0 bar of transmembrane pressure (TMP) at 1.041 m/s and 22.5 • C, in total recirculation mode to generate a quasi-steady state. The characterization process was undertaken for an operation time of 2 h to stabilize the flux through the membrane during this time. Previously each membrane was worked under compaction conditions with pure water at 5 bar for 1 h, in order to obtain a stable membrane structure. R m values were calculated using the resistance model (Equation (1)), where under this condition, there was no membrane fouling resistance (R f ) [38]: where J P is the permeate flux, µ is the viscosity of the permeate stream and R t is the total membrane resistance. R f can be understood as the result of the sum of the three main fouling mechanisms: (i) pore blockage resistance when colloids block the membrane pores, (ii) adsorption resistance as a result of foulant adsorption inside or over the membrane and, (iii) cake layer resistance as a consequence of the accumulation of particles on the membrane [35,39,40].
The evolution of permeate flux was gravimetrically measured at different time intervals according to Equation (2), where A m is the effective membrane area, m p is the total mass of permeate, ρ is the water density, and t is the filtration time.
In order to keep the feed concentration constant, both the permeate and the retained streams were continuously recirculated to the feed tank.
To evaluate the UF performance in terms of permeability the average permeate flux was calculated by the following equation [41]: where J P (t) is the permeate flux evolution over time determined by regression analysis on the experimental data, t 1 is the initial time operation (first data collected), and t M is 2 h.
To analyze the effect of the operating conditions on UF resistance, the cumulative flux decline (SFD) [41] was calculated from the following equation: where M is the number of experimental data collected and, J P (0) is the initial permeate flux measured at t 1 .
This parameter gives information about how the flux declines over the duration of the experiment (not just the difference between the initial and final permeate flux). Therefore, the higher the SFD value, the faster and more noticeable is the flux decline, indicating that the membrane fouling is more severe.
The average flux decline index (FD) provides information on the decrease of feed permeate flux throughout the experiment estimated as follows: To evaluate the UF efficiency in removing organic matter, the COD rejection rate was chosen as the response variable calculated as: where: C p and C f are the COD concentration in the permeate and feed, respectively.

Experimental Design Based on the Taguchi Method
An experimental design based on the Taguchi method was used to design the experiments. The Taguchi method applies fractional experimental designs, called orthogonal arrays (OA), to reduce the number of experiments required to determine the optimum conditions based on the results [29,30,42]. One of the important steps in the Taguchi approach is the appropriate selection of OAs, which depends on the number of control factors and their levels. The minimum number of experimental trails required in an OA is given by N min = (L − 1)·F + 1, where F and N are the number of factors and levels respectively [36,43].
As mentioned above, the four factors (parameters) chosen were the transmembrane pressure, the cross-flow velocity, the temperature and the MWCO of the membrane; and three response variables were analyzed: the average permeate flux, the COD rejection and the cumulative flux decline. The selected factors, their designated symbols and levels are presented in Table 2. For an experimental design with four factors and three levels for each factor, an L 9 (3 4 ) orthogonal array was selected. In this case, 27 runs were conducted (three repetitions at each trial condition). The design of the experiments planning matrix for the L 9 array [43] is shown in Table 3. The experiments were carried out in a randomized order. New membrane samples were preconditioned and used in the experimental runs (27 membrane samples). In this way, the lurking effect of possible irreversible fouling was avoided and the intrinsic variability of the membrane material is included in the replications. Table 3. Experimental layout using L 9 (3 4 ) orthogonal array in accordance with the Taguchi method. 1  1  1  1  1  2  1  2  2  2  3  1  3  3  3  4  2  1  2  3  5  2  2  3  1  6  2  3  1  2  7  3  1  3  2  8  3  2  1  3  9  3  3  2  1 * All experiments were carried out in a randomized run.

Experimental Trial (*) Levels
The aim of this DoE was to determine the operating parameters (factors) under which the average permeate flux and COD rejection rate achieve their maximum values, and the SFD achieved its minimum value. The Taguchi method uses a statistical measure of the process performance, called signal-to-noise ratio (S/N), which depends on the criterion for the response variable to be optimized. The S/N ratios are divided into three different categories and data sets, the larger-the better, the smaller-the-better and the nominal-the-better [16,29]. In this study, the system was optimized when the average permeate flux and COD rejection rate were as large as possible (Equation (7)), and the SFD was as small as possible (Equation (8)): The larger − the − better (S/N) = −10 log The smaller − the − better (S/N) = −10 log where Y i the response variable at each experiment and n is the total number of repetitions in a trial. The sequence of steps to be followed using the Taguchi method to optimize the UF process is shown in Figure 1. Minitab Statistical and Statgraphics Centurion XVII Software were used to analyze the Taguchi experiments and optimize the operating conditions.

Utility Concept
The implementation of the utility concept in the Taguchi method helps to obtain the best combination of operating parameters to optimize multiple response S/N ratios (MRSN) simultaneously by differentiating the relative importance (weights) of various responses [34,45,46]. In this work, it is assumed that the overall utility is the sum of the responses of the individual utilities and it can be written as [34,47]: where U(x1,x2,…xn) is the overall utility of n response variables and Ui(xi) is the utility index of i th response.
The response variables can be attributed priorities depending upon the process goals to be achieved. The priorities can be adjusted by providing a weight to the individual utility index. Therefore, by assigning weights to the response variables, the overall utility function can be expressed as: where Wi is the weight assigned to the i th response variable. It is worth noting that the assignment of weights is a purely subjective (empirical) step and depends on each experiment or process that will be carried out [48]. Therefore, in this paper, the analytic hierarchy process (AHP) method, developed by [49] was used to determine the associate weight criteria for each response variable in the multiple optimization required to calculate the overall utility index. The relative normalized weight Wi of each criterion is calculated using the AHP geometric mean method GMi on the rows in the pairwise comparison matrix, A║aij║and it can be calculated from the follow equation [49]: Minitab Statistical and Statgraphics Centurion XVII Software were used to analyze the Taguchi experiments and optimize the operating conditions.

Utility Concept
The implementation of the utility concept in the Taguchi method helps to obtain the best combination of operating parameters to optimize multiple response S/N ratios (MRSN) simultaneously by differentiating the relative importance (weights) of various responses [34,45,46]. In this work, it is assumed that the overall utility is the sum of the responses of the individual utilities and it can be written as [34,47]: where U(x 1 ,x 2 , . . . x n ) is the overall utility of n response variables and U i (x i ) is the utility index of i th response. The response variables can be attributed priorities depending upon the process goals to be achieved. The priorities can be adjusted by providing a weight to the individual utility index. Therefore, by assigning weights to the response variables, the overall utility function can be expressed as: where W i is the weight assigned to the i th response variable.
It is worth noting that the assignment of weights is a purely subjective (empirical) step and depends on each experiment or process that will be carried out [48]. Therefore, in this paper, the analytic hierarchy process (AHP) method, developed by [49] was used to determine the associate weight criteria for each response variable in the multiple optimization required to calculate the overall utility index. The relative normalized weight W i of each criterion is calculated using the AHP geometric mean method GM i on the rows in the pairwise comparison matrix, A a ij and it can be calculated from the follow equation [49]: where i,j = 1, 2, . . . , M and M is the number of factors in judgement matrix A. a ij = 1 for i = j, a ij = 1/a ij for i j.
In addition, the total sum of the weight for all the responses must be assigned to hold the following condition: For this study, as stated above, the objective was to maximize permeate flux and COD rejection rate, and minimize the SFD, simultaneously. From the utility concept, the MRSN of the overall utility value is given by Equation (13): where: W JP , W COD_rejection , and W SFD are the weights assigned to the permeate flux, COD rejection rate and SFD. It is worth mentioning that the utility function is of the higher-the-better type. If the composite measure (the overall utility) is maximized, the quality characteristics considered for the evaluation of utility will automatically be optimized (maximized or minimized) [44].

Optimal Performance Prediction
Once the optimal level of the operating conditions has been selected, it is possible to predict and verify the utility responses using the optimal parameters. The predicted response values under optimal conditions Y opt can be calculated from Equation (17) [50,51]: where m is the overall mean value of µ MRSN over nine trials; m i,j is the mean value of the quality response under optimal conditions; and p is the number of significant operating parameters that affect the UF process. The 95% confidence interval for the confirmation experiments (CI CE ) must be evaluated at the selected error level according to the following expression [43,50,51]: Membranes 2020, 10, 170 9 of 22 where F α (1,f e ) is the F-ratio at a confidence level of (1 − α) against a DOF equal to one and an error degree of freedom f e and, n eff is the effective sample size calculated as: n e f f = N 1 + (DOF of all factors used to estimate the mean) (19) where N is the number of experiments, R is the number of repetitions, and MS e is the error variance.

Analysis of Variance (ANOVA)
In order to determine the relative importance of the factors, ANOVA was employed by calculating the sum of squares (SS), degrees of freedom (DOF), mean of square (MS), associated F-test of significance (F) and percentage contribution (P%) [52].

Experimental Results
As previously mentioned, permeability tests were carried out prior to each Taguchi experiment with pure water as the feed, in order to determine the intrinsic resistance of the membrane for each membrane used, as illustrated in Figure 2.
Membranes 2020, 10, x FOR PEER REVIEW 9 of 23 where Fα (1,fe) is the F-ratio at a confidence level of (1 − α) against a DOF equal to one and an error degree of freedom fe and, neff is the effective sample size calculated as: = 1 + (DOF of all factors used to estimate the mean) (19) where N is the number of experiments, R is the number of repetitions, and MSe is the error variance.

Analysis of Variance (ANOVA)
In order to determine the relative importance of the factors, ANOVA was employed by calculating the sum of squares (SS), degrees of freedom (DOF), mean of square (MS), associated Ftest of significance (F) and percentage contribution (P%) [52].

Experimental Results
As previously mentioned, permeability tests were carried out prior to each Taguchi experiment with pure water as the feed, in order to determine the intrinsic resistance of the membrane for each membrane used, as illustrated in Figure 2. The specific resistance values obtained from the permeability test for the membranes of 10, 50 and 100 kDa were 3.46 × 10 12 , 4.56 × 10 12 and 9.88 × 10 12 m −1 , respectively.
The values for the average permeate flux, COD rejection rate and the cumulative flux decline caused by membrane fouling for each trial experiment according to the Taguchi design are shown in Table 4. The highest average permeate flux was obtained in Trial 8 ( ̅ = 95.16 L·m −2 ·h −1 ) and the lowest The specific resistance values obtained from the permeability test for the membranes of 10, 50 and 100 kDa were 3.46 × 10 12 , 4.56 × 10 12 and 9.88 × 10 12 m −1 , respectively.
The values for the average permeate flux, COD rejection rate and the cumulative flux decline caused by membrane fouling for each trial experiment according to the Taguchi design are shown in Table 4. The highest average permeate flux was obtained in Trial 8 (J P = 95.16 L·m −2 ·h −1 ) and the lowest value was obtained in Trial 1 (J P = 15.23 L −1 ·h −1 ). The corresponding average flux decline indices were respectively FD = 44.87% and a FD = 14.0%.
For some trials, significant differences on the results can be observed between replicates. The reason can be found in the inhomogeneous behavior between samples taken from the same sheet. The replication used in the experimental method aims to diminish the effect of the membrane variability on the response.  Figure 3 shows the evolution of the permeate average flux for all the trials carried out according the Taguchi design of Table 3. It can be seen that for the same MWCO, the permeate flux decreased with increasing TMP in all trials because of the membrane fouling. An increase in the TMP leads to higher accumulation of colloidal substances on the membrane surface and pore blocking. At the beginning of the process, flux declined very quickly, mainly for the trials corresponding to higher TMP and MWCO ( Figure. 4c), possibly because of the membrane pores becoming blocked more rapidly by adsorption and the accumulation of colloidal substances. Afterwards, the permeate flux continued to decline due to the growth of a cake layer on the membrane surface, until the permeate flux reached a quasi-stationary state [53][54][55][56]. For some trials, significant differences on the results can be observed between replicates. The reason can be found in the inhomogeneous behavior between samples taken from the same sheet. The replication used in the experimental method aims to diminish the effect of the membrane variability on the response. Figure 3 shows the evolution of the permeate average flux for all the trials carried out according the Taguchi design of Table 3. It can be seen that for the same MWCO, the permeate flux decreased with increasing TMP in all trials because of the membrane fouling. An increase in the TMP leads to higher accumulation of colloidal substances on the membrane surface and pore blocking. At the beginning of the process, flux declined very quickly, mainly for the trials corresponding to higher TMP and MWCO ( Figure. 4c), possibly because of the membrane pores becoming blocked more rapidly by adsorption and the accumulation of colloidal substances. Afterwards, the permeate flux continued to decline due to the growth of a cake layer on the membrane surface, until the permeate flux reached a quasi-stationary state [53][54][55][56].

Taguchi Results
The corresponding S/N ratio (in dB) calculated for the response variables at each trial are listed in Table 5.
In order to analyze the influence of each factor on the response variable, the S/N ratio for a single factor can be determined by averaging the S/N ratios at their levels. The range of the effect for each factor (Δs) is calculated as the difference between the two readings, the higher the range, the stronger

Taguchi Results
The corresponding S/N ratio (in dB) calculated for the response variables at each trial are listed in Table 5. Table 5. Signal-to-noise results (mean ± standard deviation for the three repetitions at each trial). In order to analyze the influence of each factor on the response variable, the S/N ratio for a single factor can be determined by averaging the S/N ratios at their levels. The range of the effect for each factor (∆s) is calculated as the difference between the two readings, the higher the range, the stronger the effect of the factor, in other words, it shows which parameter has the greatest effect on the response.

S/N Ratiō
The mean S/N ratio curves for each factor are shown in Figure 4. It is worth mentioning that the peak points in these plots correspond to the optimal condition.
As can be seen in Figure 4, the variations (∆s) around the mean S/N value were different for the different factors. TMP and membrane MWCO had the greatest effect on the average permeate flux as they have the steepest slope, ∆s = 7.10 and 5.95 respectively. CFV was the next one with a ∆s = 1.42, and temperature had the lowest variation around the mean S/N value, with ∆s = 0.78. In addition, from Table 5, the overall mean value was calculated as 33.59 (dB) from all the trial experiment results. It can be observed that the increase in J P was stronger when the TMP changes from 1.0 to 2.0 bar than when it changes over the range from 2.0 to 3.0 bar, this could be due to the effects of polarization and cake compaction on the membrane surface. For CFV, the slope of the line between the different levels is not the same (0.463-0.752 m/s is higher than at 0.752-1.041 m/s), but with a small variation around the J P value. Also, it can be seen for MWCO and temperature that the slopes from 10 to 100 kDa and 15 • C to 30 • C (respectively) are almost the same. Therefore, the maximum average permeate flux can be obtained for 3.0 bar, 100 kDa, 1.041 m/s CFV and high temperature (30 • C).
Under optimal COD rejection conditions, a positive and larger value of S/N is desired. In Figure 4b when comparing the S/N between different factors, it was shown that the most significant variation around the mean S/N ratios is observed for MWCO and CFV (∆s = 2.07 and 1.32 respectively). Also, it can be seen that the S/N ratio increased with TMP and CFV and decreased with MWCO. Hence maximum COD removal occurred at higher TMP and CFV (3.0 bar and 1.041 m/s), and 10 kDa. It is worth mentioning that the DCS found in the PMTE are a mixture of high and low molecular weight organic and inorganic compounds, thus the contribution of the smaller particles gives lower rejection during high MWCO UF in membranes with 50 and 100 kDa MWCO. Figure 4c shows that an increase in TMP, temperature and MWCO caused a decrease in the S/N ratio for SFD, that is to say, these factors intensified the membrane fouling effects. On the other hand, an increase in CFV induced an increase in the S/N ratio, this resulted in a decrease in the fouling effect. The highest variations around the mean S/N ratio were found for MWCO and TMP (∆s = 5.21 and 4.81). Generally, the permeate flux increased with increasing MWCO and TMP. However, under these operating conditions, DCS in PMTE can easily pass through the membrane and blocking can be observed within the pores and on the membrane surface. In addition, the highest S/N ratio for the SFD factor (−8.98 ± 0.28) was achieved in Trial 1, whereas the lowest S/N ratio (−19.43 ± 0.29) was obtained in Trial 8. The optimal conditions that minimized the SFD (lowest level of fouling) were obtained at the lowest TMP (1.0 bar), highest CFV (1.041 m/s), at temperature 15 • C and at the smallest MWCO (10 kDa). The mean S/N ratio curves for each factor are shown in Figure 4. It is worth mentioning that the peak points in these plots correspond to the optimal condition.  As can be seen in Figure 4, the variations (Δs) around the mean S/N value were different for the different factors. TMP and membrane MWCO had the greatest effect on the average permeate flux as they have the steepest slope, Δs = 7.10 and 5.95 respectively. CFV was the next one with a Δs = 1.42, and temperature had the lowest variation around the mean S/N value, with Δs = 0.78. In addition, from Table 5, the overall mean value was calculated as 33.59 (dB) from all the trial experiment results. It can be observed that the increase in ̅ was stronger when the TMP changes from 1.0 to 2.0 bar (c)

ANOVA Results
A statistical analysis of variance (ANOVA) was carried out to quantitatively determine the effect of each factor on the UF process indicators, with the aim of estimating whether the process parameters are statistically significant or not on the results responses. The ANOVA results are shown in Table 6. In order to determine the qualitative significance of each factor on the responses, Fisher's test (F-value) was employed in the ANOVA analysis. An F-value is defined as the ratio of variance due to the effect of a factor on the variance due to the inherent error in the system [57]. The F-value was compared to the critical F-value (Fcr) [52]. A calculated F-value lower than the Fcr-value means that the effect of that factor is not significant at the selected confidence level or/and it is not important in comparison with the error term. In this study, with four factors, three levels for each factor and three repetitions at each trial condition, the DOF for each factor is 2 and the DOF for the error is 18, so the Fcr-value at a confidence level of 95% is equal to 3.55. In accordance with the ANOVA table, for average permeation flux, the F-value for TMP and MWCO (114.1 and 81.52, respectively) are greater than the Fcr-value. This means that the variance of these factors is significant compared with the variance of error and they have a significant effect on the response. On the other hand, temperature and CFV had no meaningful qualitative effect on J P , as their F-values were less than the Fcr-value. Furthermore, COD rejection rate and SFD presented F-values for all factors greater than the Fcr-value which means that the effect of these factors is significant at the 95% confidence level and they have a meaningful qualitative effect on responses.
Another statistical tool that is helpful for qualitative evaluation in ANOVA is the p-value, which is used to indicate which factors had a significant effect on the responses. The smaller the p-value at an α level of significance, the more significant is the corresponding factor [58,59]. In this study, based on p-values at the 95% confidence level (α = 0.05), all the factors had a statistically significant (p-value < 0.05) effect on the COD rejection rate and SFD. For J P , CFV and temperature had a p-value higher than 0.05, thus the effect may be regarded as insignificant and it can be ignored.
The use of the percentage contribution (P%) in ANOVA analysis is helpful for the quantitative evaluation of the factorial effects of the performance indicators. The percentage contributions P% of all factors on average permeate flux, COD rejection rate, and SFD are shown in Figure 5. TMP (P% = 54.62) was the most important factor on average permeate flux, as higher pressure resulted in higher permeate flux, according to Darcy's law. MWCO (P% = 39.06) was the second most important factor, followed by temperature and CFV. For the COD rejection rate, the order of importance for the factors is as follows MWCO > CFV > temperature > TMP. In addition, MWCO and TMP (46.57% and 40.89%, respectively) were the most significant parameters on membrane fouling resistance, followed by the CFV and temperature. TMP and MWCO were the most important factors for responses. Higher TMP and MWCO resulted in higher permeate flux. However, more intensive flux decline, due to membrane fouling, occurred at higher permeate flux.
Membranes 2020, 10, x FOR PEER REVIEW 15 of 23 all factors on average permeate flux, COD rejection rate, and SFD are shown in Figure 5. TMP (P% = 54.62) was the most important factor on average permeate flux, as higher pressure resulted in higher permeate flux, according to Darcy's law. MWCO (P% = 39.06) was the second most important factor, followed by temperature and CFV. For the COD rejection rate, the order of importance for the factors is as follows MWCO > CFV > temperature > TMP. In addition, MWCO and TMP (46.57% and 40.89%, respectively) were the most significant parameters on membrane fouling resistance, followed by the CFV and temperature. TMP and MWCO were the most important factors for responses. Higher TMP and MWCO resulted in higher permeate flux. However, more intensive flux decline, due to membrane fouling, occurred at higher permeate flux. It is important to mention that the values reported due to error resulting from uncontrollable noises should be below 50% for the results to be reliable [28,29,60]. Therefore, it can be seen in Figure  5 that for average permeate flux, COD rejection rate and SFD, the error values are 4.31%, 3.07%, and 3.81%, respectively. This means that the error values for the experiment are not significant for the UF process.

Optimal Results Obtained from the Taguchi Method and Utility Concept
The aim of optimizing the process was to find the operating conditions that led to a maximum average permeate flux and COD rejection rate based on the levels that gave the highest S/N ratios for the factors (desirable values) and to minimize the SFD, that is, the levels that gave the smallest S/N ratios (adverse values).

Analysis of individual response optimization
After identifying the optimal operating conditions, the optimal responses were predicted individually using the Taguchi method and ANOVA. Table 7 shows the Taguchi prediction results for the optimal conditions for average permeate flux, COD rejection rate and SFD.
According to the Taguchi predictions, the average permeate flux at TMP 3.0 bar, CFV 0.752 m/s, at 22.5 °C and with a 100 kDa MWCO, achieves 81.20 L·m −2 ·h −1 . The COD rejection rate predicted under optimal conditions indicates a 57.92% rejection, higher than any value obtained in the DoE combinations. For the SFD under optimal conditions estimated by the Taguchi method, the minimum SFD predicted is approximately 1.80, equivalent to a ̅̅̅̅ of 8.65%. Therefore, we can see that the values of the three response variables combined are far from those values obtained experimentally (see Table 4).  It is important to mention that the values reported due to error resulting from uncontrollable noises should be below 50% for the results to be reliable [28,29,60]. Therefore, it can be seen in Figure 5 that for average permeate flux, COD rejection rate and SFD, the error values are 4.31%, 3.07%, and 3.81%, respectively. This means that the error values for the experiment are not significant for the UF process.

Optimal Results Obtained from the Taguchi Method and Utility Concept
The aim of optimizing the process was to find the operating conditions that led to a maximum average permeate flux and COD rejection rate based on the levels that gave the highest S/N ratios for the factors (desirable values) and to minimize the SFD, that is, the levels that gave the smallest S/N ratios (adverse values).

Analysis of Individual Response Optimization
After identifying the optimal operating conditions, the optimal responses were predicted individually using the Taguchi method and ANOVA. Table 7 shows the Taguchi prediction results for the optimal conditions for average permeate flux, COD rejection rate and SFD. According to the Taguchi predictions, the average permeate flux at TMP 3.0 bar, CFV 0.752 m/s, at 22.5 • C and with a 100 kDa MWCO, achieves 81.20 L·m −2 ·h −1 . The COD rejection rate predicted under optimal conditions indicates a 57.92% rejection, higher than any value obtained in the DoE combinations. For the SFD under optimal conditions estimated by the Taguchi method, the minimum SFD predicted is approximately 1.80, equivalent to a FD of 8.65%. Therefore, we can see that the values of the three response variables combined are far from those values obtained experimentally (see Table 4).

Analysis of Multi-Response Optimization
As mentioned above, in order to determine the weight for each response variable, a pairwise comparison matrix was compiled using the AHP method as presented in Table 8. Thus, the weights assigned to response variables were W JP = 0.568, W COD Rejection = 0.098 and W SFD = 0.334. The consistency ratio index (CR) is used to evaluate the consistency of AHP estimates. In this case, it was calculated as 0.021, which should be less than the allowed value of CR = 0.1, this means that the pairwise comparison matrix was considered acceptable.
The overall utility index for the µ MRSN was calculated using Equation (13) with values associated with the weights of each response, using the lager the better (S/N) and the results are presented in Table 9. ANOVA analysis was also performed for the multiple response variables using the utility concept. From Table 10 it is clear that when F-value is compared with Fcr (3.55), TMP, MWCO and CFV had a qualitatively significant effect (at a confidence level of 95%) on MRSN. The percentage contributions extracted from the ANOVA table were also used to determine the significance of each operating parameter in the process.
The P% values were arranged as follows: TMP > MWCO > CFV > Temperature. Therefore, according to the results, TMP and MWCO were the most important factors in optimizing the multi-response UF system.
The optimal operating conditions for the simultaneous response were obtained based on the criteria that both J P and COD rejection rate must be maximized and SFD should be minimized. The variation in the overall utility for the operating parameters at different levels is presented in Figure 6.   It is clear from Figure 6 that the optimal combination of operating conditions (maximum value of the overall utility) was found at the second level of transmembrane pressure (2.0 bar), the second level of cross-flow velocity (1.041 m/s), the first level of temperature (15 °C), and third level of MWCO (100 kDa).
Once the optimal levels had been selected the next step was to estimate the multi-response S/N ratio and predict the optimal values for the simultaneous optimization response, calculated using Equation (17) and presented in Table 11. Table 11. Optimal conditions for multi-response UF predicted using the utility concept.

Confirmation Experiment under Optimal Conditions
After determining the optimal operating conditions for the overall utility value and the significance of factors, validation experiments (for multi-responses) were carried out at the optimal levels in order to validate the predicted UF responses suggested using the Taguchi method with utility concept [34,51].
The observed permeate flux results as a function of time under optimized conditions during the UF of PMTE are plotted in Figure 7. As described previously, the flux decline was mainly the result of two phenomena, pore blocking and cake layer formation, which mostly occurred in the first hour It is clear from Figure 6 that the optimal combination of operating conditions (maximum value of the overall utility) was found at the second level of transmembrane pressure (2.0 bar), the second level of cross-flow velocity (1.041 m/s), the first level of temperature (15 • C), and third level of MWCO (100 kDa).
Once the optimal levels had been selected the next step was to estimate the multi-response S/N ratio and predict the optimal values for the simultaneous optimization response, calculated using Equation (17) and presented in Table 11.

Confirmation Experiment under Optimal Conditions
After determining the optimal operating conditions for the overall utility value and the significance of factors, validation experiments (for multi-responses) were carried out at the optimal levels in order to validate the predicted UF responses suggested using the Taguchi method with utility concept [34,51].
The observed permeate flux results as a function of time under optimized conditions during the UF of PMTE are plotted in Figure 7. As described previously, the flux decline was mainly the result of two phenomena, pore blocking and cake layer formation, which mostly occurred in the first hour of the process [21,29]. During the first 30 min of the UF, the flux decreased by 22.63%. Furthermore, immediately after pore blockage, the permeate flux continued to decline due to the formation and growth of a cake layer until the system approached the quasi-steady state. At the end of the process (after 2 h), the final permeate flux was about 67.0 L·m −2 ·h −1 and flux decline was around 39.72%, which confirms that the membrane fouling took place with a higher rate in the first 30 min and at a slower rate when the system had achieved a steady state. Therefore, the observed experimental values of average permeate flux and cumulative flux decline were about 81.15 L·m −2 ·h −1 and 6.01 (SFD equivalent to a FD value of 28.96%).
Membranes 2020, 10, x FOR PEER REVIEW 18 of 23 of the process [21,29]. During the first 30 min of the UF, the flux decreased by 22.63%. Furthermore, immediately after pore blockage, the permeate flux continued to decline due to the formation and growth of a cake layer until the system approached the quasi-steady state. At the end of the process (after 2 h), the final permeate flux was about 67.0 L·m −2 ·h −1 and flux decline was around 39.72%, which confirms that the membrane fouling took place with a higher rate in the first 30 min and at a slower rate when the system had achieved a steady state. Therefore, the observed experimental values of average permeate flux and cumulative flux decline were about 81.15 L·m −2 ·h −1 and 6.01 (SFD equivalent to a ̅̅̅̅ value of 28.96%). The total resistance ( ) at the end of the 2-h experiment under optimal conditions was 1.13 × 10 13 m −1 , which is sum of the intrinsic membrane resistance (Rm = 3.40 × 10 12 m −1 ) and fouling resistance (Rf = 7.91 × 10 12 m −1 ).
In addition, the membrane surface morphologies were observed by FESEM. Figure 8 shows the images of the membrane (PES 100 kDa) before and after the UF experiments were carried out. As can clearly be seen, before UF there is no blocking on the pores and no cake layer on the membrane surface. Figure 8b,c shows the surface of the membrane fouling after 30 min and after 2 h. In both cases, the images show the existence of pore blocking due to DCS adsorption within the membrane pores and sediments deposited on the surface (cake layer) acting to resist the UF [61]. Figure 8(d) shows the morphologies of the fouling sediments on the membrane. Furthermore, it can be seen that the membrane was indeed fouled after 30 min filtration. However, the FESEM images of the membranes after filtration, at 2 h, were highly similar to the membranes after 30 min filtration. Therefore, it may be concluded that the permeate flux decline might result from the pore blocking as opposed to the formation and growth of the cake layer on the membrane.
Additionally, to verify the permeate applicability for paper mill reuse, the physical and chemical properties of the treated effluent obtained under optimal operational conditions (PES 100 kDa membrane at TMP = 2.0 bar, CFV = 0.752 m/s, and T = 15 °C) were compared with treated paper mil effluent used in this study as feed solution ( Table 1). The results obtained for the physical-chemical parameters are given in Table 12. From the results obtained, all properties showed high retention efficiencies and proved the effectiveness of the UF under optimal conditions. The total resistance (R t ) at the end of the 2-h experiment under optimal conditions was 1.13 × 10 13 m −1 , which is sum of the intrinsic membrane resistance (R m = 3.40 × 10 12 m −1 ) and fouling resistance (R f = 7.91 × 10 12 m −1 ).
In addition, the membrane surface morphologies were observed by FESEM. Figure 8 shows the images of the membrane (PES 100 kDa) before and after the UF experiments were carried out. As can clearly be seen, before UF there is no blocking on the pores and no cake layer on the membrane surface. Figure 8b,c shows the surface of the membrane fouling after 30 min and after 2 h. In both cases, the images show the existence of pore blocking due to DCS adsorption within the membrane pores and sediments deposited on the surface (cake layer) acting to resist the UF [61]. Figure 8d shows the morphologies of the fouling sediments on the membrane. Furthermore, it can be seen that the membrane was indeed fouled after 30 min filtration. However, the FESEM images of the membranes after filtration, at 2 h, were highly similar to the membranes after 30 min filtration. Therefore, it may be concluded that the permeate flux decline might result from the pore blocking as opposed to the formation and growth of the cake layer on the membrane.
Additionally, to verify the permeate applicability for paper mill reuse, the physical and chemical properties of the treated effluent obtained under optimal operational conditions (PES 100 kDa membrane at TMP = 2.0 bar, CFV = 0.752 m/s, and T = 15 • C) were compared with treated paper mil effluent used in this study as feed solution ( Table 1). The results obtained for the physical-chemical parameters are given in Table 12. From the results obtained, all properties showed high retention efficiencies and proved the effectiveness of the UF under optimal conditions. Membranes 2020, 10, x FOR PEER REVIEW 19 of 23 (a) (b).
(c) (d).  The optimum predicted results at the 95% confidence interval calculated using Equation (18) and the observed experimental results for the response variables are given in Table 13.
The observed multi-response of the overall utility falls within the 95% confidence interval for the optimal range of the response variables. In addition, it is clearly observed that the deviation between predicted and experimental results is very small, which confirms that the Taguchi method and utility concept can be used to predict the multi-response UF for any parametric combination, while individual optimization don't got good predictions.  The optimum predicted results at the 95% confidence interval calculated using Equation (18) and the observed experimental results for the response variables are given in Table 13.
The observed multi-response of the overall utility falls within the 95% confidence interval for the optimal range of the response variables. In addition, it is clearly observed that the deviation between predicted and experimental results is very small, which confirms that the Taguchi method and utility concept can be used to predict the multi-response UF for any parametric combination, while individual optimization don't got good predictions.

Conclusions
In this study, the Taguchi method, utility concept and ANOVA analysis were used as statistical tools to investigate the effects and significance of four operating parameters and to optimize the UF process with respect to average permeate flux, COD rejection rate and cumulative flux decline.
ANOVA was used to determine the most significant factors affecting the response variables. From the percentage contribution, the order of importance of each factor in maximum J P was TMP > MWCO > T > CFV; for maximum COD rejection rate it was MWCO > CFV > T > TMP; and to achieve the minimum SFD: MWCO > TMP > CFV > T.
The optimal UF operating parameters, based on the Taguchi method and utility concept, were found at TMP (2.0 bar), CFV (1.041 m/s), temperature (15 • C) and MWCO (100 kDa). Under these optimal conditions, J P , COD rejection rate and SFD resistance of 81.15 L·m −2 ·h −1 , 43.90% and 6.01 (around and FD value of 28.96 %), respectively, were obtained and they were within of the predicted range at the 95% confidence interval.
Measurements of turbidity, COD and particle size in the permeate showed a significant decrease 3.21 to 0.0002 NTU, 146 mg/L to 81.8 mg/L and 458-1281 nm to 12.71-24.36 nm, respectively, which confirms a substantial reduction in colloidal compounds. Therefore, it can be said that UF is suitable for removing dissolved and colloidal substances from wastewater effluents from recycled paperboard manufacturing.
Finally, we can say that the Taguchi method and utility allow membrane conditions for the P & P wastewater treatment to be optimized using a reduced number of experiments. The methodology used in this study could be used as a guideline for operating UF systems applied as a tertiary treatment for paperboard mill treated effluents under optimal conditions.