Maximizing Adsorption Involving Three Solutes on Enhanced Adsorbents Using the Mixture-Process Variable Design.

Unmodified (UN), acid-treated (AT) and microwave-acid-treated (MAT) activated carbons were optimized for their solute removal efficacies by adjusting feed mixture compositions and process conditions. Acetaminophen, benzotriazole, and caffeine were used either individually or as binary/ternary mixtures in this study. The process conditions considered were the pH, adsorbent dosage, and type of adsorbent. Experimental responses such as total adsorbent loading (q total) and total percentage removal (PRtotal) were fitted with empirical models that had high adjusted R 2 (>0.95), insignificant lack of fit (p-value > 0.22), and high model predictive R 2 (>0.93). Mixture compositions of the feed were found to interact significantly not only among themselves but with process variables as well. Hence, adsorption optimization must simultaneously consider mixture as well as process variables. The conventional response surface methodology for mixtures, termed as ridge analysis, optimizes mixture compositions at specified values of process variables. An improved steepest ascent method which considers mixture and process variables simultaneously was developed in this work. This could track the path of steepest ascent toward globally optimal settings, from any arbitrary starting point within the design space. For the chosen adsorbent, optimal settings for feed mixture compositions and pH were found to change along this steepest ascent path. The feed compositions, pH, and adsorbent dosage identified for maximum adsorbent utilization were usually quite different from those identified for maximum total percentage removal. When both these objectives were optimized together, the most favorable compromise solutions for q total and PRtotal were, respectively, 264.1 mg/g and 43.4% for UN, 294.9 mg/g and 52.5% for AT, and 336.6 mg/g and 55.9% for MAT.

Introduction

Adsorption is a popular separation-purification process that removes different solutes from the fluid phase to various extents using a solid adsorbent. Even though simultaneous adsorption of multiple solutes is practically relevant, such studies are relatively less when compared to single-component adsorption in the literature. In multicomponent systems, the affinities of different solutes toward the adsorbent may be quite different from those of individual solutes.[1,2] A recent study by Chen et al. (2022)[3] indicated that when chromium and humic acid were present together in wastewater, synergetic interactions led to both being adsorbed to a greater extent by the powdered activated carbon (AC) adsorbent than when they were present alone. Onaga Medina et al. (2021)[4] studied the binary adsorption of diclofenac and caffeine on AC. They observed that the presence of either solute in the mixture favorably influenced the adsorption of the other.

The presence of multiple solutes are typically encountered in wastewater containing multifarious pollutants of different chemical structures and household gray water that typically has chemical compounds sourced from washing detergents, personal care products, medicines, beverages, and dish wash soaps.[5] Thus, understanding and optimizing multicomponent adsorption is vital to both industry and society.

AC is one of the preferred adsorbents for wastewater treatment[6] owing to its versatility, low cost, and ability to be made from numerous locally available organic sources.[7] AC is a complex material whose physicochemical properties may be altered to modify its adsorption capacities and affinities for different compounds.[8] Alterations may be carried out through thermal treatments that include conventional or microwave heating and chemical treatments using oxidants, acid, base, or other reagents.[9,10] Conventionally, sulfuric acid treatment has been proven to introduce oxygen-containing functional groups over the surface of AC.[11,12] Furthermore, microwave treatment may be coupled with chemical modifications either simultaneously[13,14] or sequentially.[15,16] Hence, there is considerable scope for further enhancing adsorbent performances through different post activation treatments. These enhanced adsorbents may be made to perform optimally by identifying suitable feed mixture compositions and/or operating conditions that are conducive to adsorption as there may not be universal conditions for different adsorbent types as well as different solutes.

Solutes’ feed compositions in the mixture, nature of the adsorbent, and operating conditions such as pH, temperature, and adsorbent dosage influence the loading of solutes on the adsorbent. The compositions of solutes in the aqueous solution are termed as mixture variables, while the nature of the adsorbent, pH, and adsorbent dosage are termed as process variables.

Design of experiments (DOE) provides valuable insights while economically varying the levels of variables or factors.[17] A few studies dealing with multicomponent mixtures have varied the total solute concentrations at two levels and studied their role on the adsorbent’s loading capacity.[18,19] A more sophisticated DOE strategy entails a mixture design which investigates the influence of different relative proportions of the constituents while keeping the total mixture concentration constant.[20,21] Factorial design and mixture design approaches are compared in Table S1 in the Supporting Information.

Combining the mixture design with process variables results in a more inclusive mixture-process variable (MPV) design. This approach has been considered by few studies to investigate the simultaneous effect of process variables and solute composition on the adsorption process.[22,23] MPV models that correlate responses with the factors and their interactions have to be validated against new experimental data. In the literature, the validated models have been optimized either for single response[24,25] or for multiple objectives.[26]

Zolgharnein et al.[22,23] utilized the MPV approach to find the optimal initial concentration, pH, and adsorbent dosage that resulted in maximum biosorption of three heavy metals. The MPV approach is becoming increasingly popular and has been applied in other fields as well, for example, in food processing. Nasehi et al.[27] optimized the formulation of spaghetti using the MPV design that consisted of three mixture variables and two extrusion-related process variables. They made inferences from the mixture surface and contour plots on a large number of nutritional and sensory properties. Kashaninejad et al.[28] applied the MPV design with two mixture variables and one process variable to optimize the production of labane. Their analysis used the desirability optimization approach using Design Expert (Stat-Ease, Inc., Minneapolis). These studies indicate that the inclusion and quantification of interactions in the MPV model facilitate reliable model development, optimization, and deeper insights into the process.

However, interpreting the effect of process conditions as well as mixture compositions on adsorbent capacity[19,23] is not straightforward. Complications arise when different variables affect the overall adsorption either individually or through their interactions with one another in different possible combinations. These include not only interactions among compositions and process variables separately but also those between them as well. Often these interactions may exert even more influence than the individual variables,[20,29] but such combined MPV interaction studies in multicomponent adsorption are scarce in the literature.

Ridge analysis procedures have been detailed for the steepest ascent toward global optimum involving only mixture variables[30] or only process variables.[31] These paths represent the loci of locally optimal solutions, which may be potentially considered if the global optimum conditions are not feasible to the operating treatment facility. However, to our knowledge, tracing the path of the steepest ascent from any arbitrary starting combination of MPVs toward the globally optimal adsorption performance has not been detailed for multicomponent adsorption problems.

Three model compounds that considerably vary in their properties and are representative of their respective groups have been chosen for this study. These are acetaminophen (or paracetamol), benzotriazole (a chemical used in washing machine detergents), and caffeine (an important beverage stimulant). They differ in their distribution coefficient values (log Kow) and dissociation constants (pKa), which are listed along with other physicochemical properties as Table S2 in Supporting Information. These solutes also had high frequency of detection in the influents of WWTP as well as high persistence in the environment.[32,33]

Based on the above, the objectives of the present multicomponent adsorption study are as follows:

Using a compact MPV design, investigate the main and interaction effects of mixture composition, pH, adsorbent dosage, and adsorbent type on the total adsorbent loading qtotal of the solutes

Identify the loci of locally optimal maximum total adsorbent loadings and percentage removals (PRtotal), individually as well as their weighted sum objective, en route to their respective global maximum values from any initial composition and process condition. The process and composition variables at each of the locus points are also to be identified when tracking the steepest ascent path toward the global maximum.

The novel aspects in the present study are given below:

Maximizing (a) total solute uptake by the adsorbent qtotal and (b) percentage removal PRtotal in multicomponent adsorption considering mixture compositions and process variables simultaneously.

Synergetic and antagonistic interaction effects between mixture variables, between process variables, and also between mixture and process variables were quantified.

The theory-based ridge analysis optimization procedure for mixtures[29] was considerably improved in our work for handling both mixture and process variables simultaneously, and this represents a valuable new contribution.

Globally optimal conditions were separately identified for maximum adsorbent loading qtotal and maximum percentage removal PRtotal. In addition, locally maximum conditions were also identified along the path of steepest ascent toward the global maximum in our improved method. This also gives the flexibility to carry out the adsorption in a locally optimal manner at specified compositions of the feed as dictated by field conditions. The best AC, pH, dosage, and mixture feed compositions are known at each local optimum solution.

Materials and Methods

Materials

Commercial AC procured from Active Char Products Pvt. Ltd. Edyar, Kerala, was washed, dried in a vacuum oven for 24 h at 110 °C, sieved to 0.425–0.5 mm, and labeled as UN. The three model solutes used, acetaminophen (ACT), benzotriazole (BTA), and caffeine (CAF), were of analytical grade purity and were bought from SD Fine-Chem Ltd., Loba Chemie Pvt. Ltd., and HiMedia Laboratories, respectively. Concentrated sulfuric acid (98% w/w) procured from Sisco Research Laboratories Pvt. Ltd., Mumbai, was of analytical grade and was diluted to 1 M using ultrapure water from the purifier of Evoqua Water Technologies, Pennsylvania, US. Acetonitrile of HPLC grade purchased from Finar Pvt. Ltd., Mumbai, was used to prepare the HPLC mobile phase. The aqueous solutions were prepared using ultrapure water.

Sulfuric Acid and Microwave Modification of AC

As per Li et al.,[34] 20 g of corresponding AC (UN) was stirred in 400 mL of 1 M H2SO4 at 400–500 rpm for 3 h in a constant water bath at 60 °C. This treated carbon was washed in distilled water until pH increased to that of the washing media. The washed carbon was dried at 110 °C for 24 h in a vacuum oven. This carbon was labeled as acid-treated (AT). For microwave treatment, a T-neck containing cylindrical quartz tube (ID, 2.5 cm; height, 30 cm) was inserted inside a microwave oven (MW73AD, Samsung) from the top. 10 g of UN carbon was added to the cylindrical quartz tube that was purged continuously with N2 and exposed to 450 W microwave power for 20 min. The optimal microwave power and exposure time were determined from preliminary studies as 450 W and 20 min, respectively. This microwave-treated carbon was further subjected to 3 h of acid treatment as mentioned above and labeled as microwave-acid-treated (MAT) carbon.

Carbon Characterization

The ACs (UN, AT, and MAT) were characterized using BET, pHpzc, and FTIR studies. BET isotherm experiments were carried out using nitrogen as analysis gas by Sprint Testing Solutions, Mumbai, using a Quantachrome ASiQwin instrument after degassing the sample at 200 °C for 10 h. To acquire the Fourier transform infrared (FTIR) transmission spectrum, the AC sample was crushed, mixed with KBr, and pelletized. This pellet was scanned using a PerkinElmer FTIR spectrometer in the range of 450–4000 cm–1 with 1 cm–1 increments. The FTIR spectra were determined at Sophisticated Analytical Instruments Facility (SAIF), IIT, Madras. Scanning electron microscopy (SEM) images were obtained from a Hitachi S 4800 after gold-sputtering the carbons. The point of zero charge (pHpzc) for different adsorbents was estimated using the salt addition method described by Gil et al.[35]

Batch Adsorption Procedure

The total initial concentration was fixed at 700 mg/L in the experimental design. Subject to this constraint, 100 mL solutions containing different proportions of the three solutes were prepared. The initial pH, measured using a Eutech pH 700 pH meter, was adjusted using 0.1 N HCl and 0.1 N NaOH. Appropriate carbons (UN, AT, and MAT) were added to the samples and shaken at 130 rpm and 27 °C for 72 h in an orbital shaker, after which the equilibrium concentration was detected using HPLC.

The concentration of the three solutes in liquid was measured simultaneously by an isocratic HPLC procedure using a C18 column (KyaTech Japan) attached to the Jasco 2010 equipment with a photodiode array detector. The mobile phase was 0.01 M KH2PO4 at pH 3 (80% v/v) mixed with acetonitrile (20% v/v) flowing at 0.8 mL per minute. The absorbance was observed only at 260 nm after confirming its efficacy by comparing the values at 243, 260, and 273 nm during preliminary studies. The mobile phase was mixed well, vacuum-filtered, ultra-sonicated, and cooled to room temperature before being pumped across the C18 column. The experimental errors in total adsorbent loading and total percentage removal were estimated to be ±1.9 and ±1.5%, respectively.

MPV Design and Model Equation

Design Expert 11 (Stat-Ease, Inc. Minneapolis) software was used to carry out the experimental design, model development, and subsequent analysis. Optimization was carried out using MATLAB 2018b (The MathWorks, Inc. Natick, Massachusetts), once the MPV model was validated.

The composition space was represented in the form of an equilateral triangle, and single, binary, and ternary mixtures were represented in it. The sum of the three concentrations was constrained to be 700 mg/L (eq )

In order to investigate this region completely, vertices, centers and thirds of edges, axial check blends, interior check blends, and the overall centroid were chosen as the 19 candidate points. These points are depicted in Figure S1 of Supporting Information. The distance-based optimality criterion was preferred to disperse the selected points uniformly throughout the design space. Using this approach, design points that may lead to unusual combination of the factors such as low non-zero concentrations could be avoided.[36]

In addition, three process variables, namely, (a) pH at three discrete levels (3, 6.5, and 10), (b) adsorbent dosage as a continuous variable assigned between 0.6 and 1.2 g/L, and (c) carbon type as a category variable with three levels (UN, AT and MAT), were chosen. For the suggested 91 runs, solute concentrations in the liquid were measured after 72 h of equilibration, using HPLC, and were used to calculate different responses. The response values along with actual and coded values of the factors are provided as Table S3 in Supporting Information.

The coefficients of an empirical model that relates the selected response to a combination of the input variables may be estimated by linear regression. Significant terms in the quadratic by the quadratic model (eq 2) were identified from the contribution of each term to the regression sum of squares using the analysis of variance (ANOVA) technique. The overall response Roverall may be defined as a product of two terms, that is, RMRP. RM and RP are defined in eqs and 2b, respectively.

Therefore,

The variables A, B, and C indicate the initial concentrations of ACT, BTA, and CAF solutes (in mg/L), respectively, while D, E, and F indicate the pH, dosage (mg/100 mL), and type of carbon, respectively. It may be seen that the order of both the composition and process variable models is 2. The empirical model with only significant and necessary terms was used in further analysis. The most suitable form of eq 2 was identified based on criteria such as high-adjusted R2, high-predicted R2, insignificant lack of fit, and the linear normal probability plot of the residuals. The model selection option in Design Expert 11 was used for this purpose.

Optimization of the Response Model

The MPV models were validated experimentally at newly chosen random conditions and at optimal conditions suggested by Design Expert. For validation, pH was considered to be a discrete variable with values of 3, 6.5, and 10. For the validated MPV models, global optima of qTotal and PRtotal for the three adsorbents were identified using the particle swarm optimization (PSO) routine of MATLAB R2018b. Here, pH was considered as a continuous variable with a range of 3–10.

The empirical model (eq 2) with process variables specified beforehand may also be optimized by a response surface methodology termed ridge analysis. It is the method of steepest ascent toward the optimum solution in the composition space and can be applied only for second-order models.[31] Owing to its inherent limitations, novel constrained optimization approaches were also developed.

Ridge Analysis for Mixture Designs

The evolution of the maximum response and the corresponding optimal settings of factors may be captured by ridge analysis in a single graphical plot.[30] The ridge analysis can be carried out only for quadratic equations. The complete theory and equations for the conventional ridge analysis[29] are given in the Supporting Information section.

When the process conditions such as pH (D), dosage (E), and type of carbon (F) are specified beforehand, the model (eq 2) reduces to a quadratic form with only the mixture variables as below.

In ridge analysis, a particular focus f is chosen in the ternary composition space as a starting point. Concentric circles of increasing radii R are constructed with the focus f as the center (Figure a). The ridge analysis procedure enables the identification of an optimal response, which is constrained to lie on each circle centered on the focus f. This focus may be fixed at the centroid of the triangle or at any arbitrary point on any one of the three binary edges or even within the triangle (Figure a). Also plotted qualitatively in this diagram is the locus of points, where the response is the maximum on each circle. The limitation of this method is that the process variables have to be specified a priori and only the mixture variables are allowed to vary during the optimization exercise. The ridge analysis has to terminate once the optimal compositions are identified beyond the feasible composition space, which is the triangular domain including its boundaries (Figure b).

Figure 1

(a) Circles originating from the focus arbitrarily located on the BTA-CAF edge is considered to illustrate the ridge analysis. The optimum may eventually lie outside the triangular composition space with increasing distance from the focus f. (b) Truncated circles are only considered after imposing the composition constraint that the search region should not lie outside the composition space.

Optima predicted outside the design space are not reliable as these may involve non-realizable compositions such as negative concentrations. Even if physically meaningful, the optimum predicted may not be reliable as the model begins losing its predictive capability rapidly, when the compositions lie outside the range of values used to develop it. Further, this analysis assumes prior specification of process variables. Hence, as a next step to increase the utility of ridge analysis, it is required to respect the constraints imposed by the composition bounds and simultaneously consider the process variables along with the mixture variables in the optimization exercise.

Circular Constrained Optimization for Mixtures and MPV Designs

To explore the ternary composition space completely including its borders, regardless of where the focus is located, a new constrained numerical optimization procedure is proposed. The objective function equation and associated constraints of the new circular constrained optimization procedure are given below in eq 4.

Equation is the objective function that involves maximization of qtotal.

Equation is the mixture constraint which stipulates that the total concentration (be it a single, binary, or ternary component) of the mixture should not exceed 700 mg/L.

Equation refers to the composition constraint which fixes the range of each solute composition from 0 to 700 mg/L (both inclusive).

Equation fixes the range of solution pH from 3 to 10 (both inclusive).

Equation fixes the range of adsorbent dosage from 0.6 to 1.2 g/L (both inclusive).

Equation is the nonlinear equality constraint, which is termed as the circular constraint.

With constraints from eqs –4e and the circular constraint given by eq , it is ensured that for any given radius of the circle, the search domain for maximum qtotal lies within the triangle or on the triangle’s edges (Figure a).

As discussed in Section , the ridge analysis fails once the optimum value on the circle lies outside the ternary composition domain. However, in such cases, it is necessary to identify the next best possible optimum response which lies within or on the borders (i.e., triangular edges) of the composition space. Based on this requirement, the individual composition constraints (eq ), as shown above, were introduced in the new optimization scheme termed as the circular constrained optimization method. The conventional ridge analysis inherently utilized only the circular constraint (eq ) and the total mixture composition constraint (eq ) but not the individual composition constraints (eq ). Further details on the ridge analysis may be seen in the Supporting Information section.

In the circular constrained optimization method, either mixture variables alone (after setting process variables to fixed values) or both mixture and process variables can be varied. The latter option is not possible in ridge analysis. The second difference is that composition constraints (eq ) in the circular constrained method limit the composition space to within or the edges of the triangular composition domain. Hence, negative compositions or out of the triangular domain solutions, as encountered in conventional ridge analysis, can be ruled out. The circular constrained method is shown in Figure b.

The circular constrained optimization strategy was carried out on both mixture (eq ) and MPV models (eq 2). Evolution of the maximum response could be traced with increasing radial distances from the focal point for both these designs.

To ensure that the path of steepest ascent in MPV optimization passed through the globally optimal response, the circular constrained optimization method as discussed above had to be improved. A new strategy termed cyclic optimization was used. Rather than simultaneously optimizing both process and mixture variables, they were optimized separately in two stages. The algorithm is described below.

Select the type of adsorbent.

Choose any location in the triangular composition space or on its edges as the starting point. This is called the focus f.

Calculate the distance Rfinal between the global optimal (termination) point and the focus.

Start at the focus.

For the chosen adsorbent, provide initial guesses for the following:

Feed concentrations of ACT, BTA, and CAF solutes.

Process variables (pH and dosage).

These initial guesses are used only in the first iteration.

If the distance from focus is ≤Rfinal, go to step G, else go to step M.

Define the two stages. Let the mixture variable optimization section be referred to as stage 1 and the process variable optimization section as stage 2.

In stage 1, maintain the process variables at previous iteration’s stage 2 values. Substitute them in eq . Implement constrained optimization in MATLAB (using eqs , 4c, and 4f) and find optimal compositions for only mixture variables.

Send the mixture variable values and process variables values from stage 1 to stage 2.

In stage 2, keep the mixture variables’ values from current iteration’s stage 1 constant and use in eq . Implement constrained optimization in MATLAB (using eqs and 4e) and find optimal settings for process variables.

Compare the difference between the corresponding values of all the variables stored in the two stages for that iteration. Calculate the Euclidean distance between the solutions in the two stages. If this value is not below the specified low tolerance value, go back to step G, else go to step L.

Increase the distance from the focus by a small value. Go to step F.

Steepest ascent path has terminated at the global maximum response value. End iterations.

For the specified adsorbent, this strategy results in two stages. The first stage optimizes three composition variables, while the second stage optimizes two process variables as explained above. This considerably eases the optimization search and facilitates the identification of correct local optimum solutions on the concentric circles that are centered at the focus. The steepest ascent path now includes the global optimum response irrespective of the location of the focus f and satisfies the feasible MPV domain.

Results and Discussion

Analysis of MPV Design

In a preliminary analysis, the experimental results (given in Table S3) were evaluated by considering only one factor (also called the variable) rather than all factors simultaneously. The penalty for this simplification was the spread in experimental data at each level of factor that is considered in isolation. Obviously, the factors may mutually interact and influence each other in affecting the adsorption, and this contribution is analyzed in detail in subsequent sections. However, this one variable analysis revealed useful initial results.

Effect of Variation of Each Factor Assuming the Absence of Other Factors

In the present study, the total solute loading on the adsorbent (qtotal) was obtained by summing the individual adsorbent loading, which in turn is defined asHere, ma is the mass of the adsorbent and VL is the liquid volume. The ratio ma/VL is termed here as dosage. The experimental qtotal is plotted as a box plot for each process variable, namely, AC type, pH, and adsorbent dosage in Figure .

Figure 2

Experimental responses of total loading (qtotal mg/g) for each process variable revealing the uniqueness of AC type and pH = 10. The horizontal line within each box represents the median adsorbent loading.

In a preliminary analysis presented in Figure , we present the results of varying the setting of a single process variable at a time such as (a) varying only the adsorbent type, (b) varying only the pH, and (c) varying only the adsorbent dosage. At each setting of type of adsorbent in Figure a, there is a wide variability in responses shown by the box plots as their whiskers and even outliers. This variability indicates that other variables and their interactions are also responsible for influencing the qtotal responses. Even with the variability in responses, this preliminary single variable analysis shows strong influence of the type of adsorbent and pH. From the median values (central line within the boxplot), it may be observed that pH (esp. 6.5 and 10) and type of adsorbent have a significant influence on qtotal, while adsorbent dosage does not have such a strong effect. 95% confidence intervals for difference in the means between the two settings of the variable considered are plotted in Figure S2. If the 95% confidence interval encompasses 0, then the difference between the responses at the settings compared is statistically insignificant.

Table summarizes the BET surface area, pHpzc, and FTIR characterizations of the ACs utilized in the present study. The FTIR spectra of the three carbons are presented in Figure S3a, and the characteristic peaks are discussed in the Supporting Information. The scanning electron microscopy (SEM) images of the three ACs are shown in Figure S3b. Well-developed pore structures in AT and MAT ACs may be seen in the images. The increased adsorption by AT carbon when compared to UN (Figure a) could be due to the increased BET surface area and pore volume. The increased adsorption by MAT could be due to its pHpzc value (Figure a and Table ). This indicates that favorable chemistry in MAT may outweigh the enhancement in total adsorption capacity due to the increased surface area in AT carbon.

Table 1

Summary of the Physicochemical Characteristics of the Three Activated Carbons

carbon	UN	AT	MAT
BET surface area (m2/g)	837.2	996.2	865.3
pore volume (cm3/g)	0.281	1.286	0.297
pHpzc	7.2	5.5	6.3
surface charge at pH 3	+	+	+
surface charge at pH 6.5	+	–	neutral
surface charge at pH 10	–	–	–
FTIR	oxygen-containing functional groups (C=O, O–H, −COOH) are present

Along these lines, Galhetas et al.[37] have commented on the role of pore structure and surface chemistry in adsorption. For the solutes, they observed that pore dimensions and surface chemistry determine the affinity of ACT and CAF, respectively, toward AC.

The pKa of BTA is 8.2, hence it completely dissociates into negatively charged species at pH 10. For ACT and CAF, as their pKa values are 9.5 and 10.4, respectively, their neutral species dominate the aqueous solution at pH 10. Based on the pHpzc plots, the net surface charge of carbons was found to be negative at pH 10. Thus, the low adsorption of BTA at pH 10 is due to the electrostatic repulsion between this solute and the negatively charged carbon surface.

Analysis of Variance

ANOVA was first used to identify the significant factors and the interactions before finalizing the regression equation.[38,39] ANOVA provided as Table identified the statistically significant terms in model eq . As indicated in Table of the paper, the pure error variance (mean square error) was about 91, while the mean square of most effects were in the range of 400–24000. This indicates that the variation from these effects were much higher than the variation caused by random errors. The random errors were estimated from experiment replicates and are summarized as pure error in Table . The main factors and interactions that had p-values[29] below 0.05 were considered to be statistically significant and were included in the final model. The most suitable form of the empirical regression model that led to best statistical parameters was chosen as mentioned in Section . The unequal effects of the factors and their interactions were shown by the differences in their associated p-values.

Table 2

Analysis of Variance of qtotal Where A (ACT), B (BTA), and C (CAF) Are Mixture Concentrations, D Is the pH, E Is the Adsorbent Dosage, and F Is the Type of Carbon

source	sum of squares	DFa	mean square	F-value	p-value
model	1.621 × 105	26	6236.41	71.31	<0.0001	significantb
BD	24486.42	1	24486.42	280.01	<0.0001
CF	18277.57	2	9138.78	104.50	<0.0001
AF	11137.46	2	5568.73	63.68	<0.0001
BF	10566.43	2	5283.21	60.41	<0.0001
BC	8348.38	1	8348.38	95.46	<0.0001
BD2	7775.53	1	7775.53	88.91	<0.0001
ABD	1694.66	1	1694.66	19.38	<0.0001
AC	1634.76	1	1634.76	18.69	<0.0001
BDF	1823.21	2	911.60	10.42	0.0001
CEF	1546.80	2	773.40	8.84	0.0004
BE	1051.66	1	1051.66	12.03	0.0009
ADF	1349.50	2	674.75	7.72	0.0010
AB	720.29	1	720.29	8.24	0.0056
AD	678.57	1	678.57	7.76	0.0070
linear mixture	873.96	2	436.98	5.00	0.0096
CDE	492.75	1	492.75	5.63	0.0206
BCD	424.74	1	424.74	4.86	0.0311
CD	165.47	1	165.47	1.89	0.1737	insignificant
lack of fit	4594.28	53	86.68	0.9511	0.5830	insignificant
CE	9.12	1	9.12	0.1043	0.7478	insignificant
residual	5596.78	64	87.45
pure error	1002.50	11	91.14
total	1.677 × 105	90

DF is degrees of freedom.

The bold-faced terms are statistically significant (p-value < 0.05).

When building the model, not all insignificant terms were removed. Insignificant terms whose higher order combinations are significant were retained to maintain the model hierarchy. For instance, in Table , we may observe that CE and CD terms are insignificant, yet their combination, CDE, is significant (p-value = 0.0206). Hence, either all these three terms must be removed or retained together to maintain model hierarchy.

Removing these three terms led to lower adjusted R2 and hence the model was not trimmed down further. Since the lack of fit is insignificant, it is not necessary to add more terms. Analyzing the p-values in Table , we observed that the linear mixture effects representing contributions from individual solute concentrations were significant (p-value < 0.01). However, these contributions were relatively weaker when compared to the binary interaction terms, that is, combinations of two factors as the latter had p-values smaller than 0.0001.

The interactions of different solutes with the type of carbon (e.g., AF or BF or CF) and pH (e.g., BD) populate the highly significant binary terms (p-values < 10–4) and contribute considerably to the variability in the process response. A strong conclusion that may be made is that the type of carbon shows significant interactions with all three solutes, suggesting that the modification procedures affect the adsorption of individual solutes in different ways. The BTA concentration (B) strongly interacts with pH (D) as the quadratic dependence on pH (BD2) manifests only when BTA is present. The solutes interact significantly among themselves as well, as the terms AB, BC, and AC are significant. However, their nature of interaction may be further subject to pH changes since terms ABD and BCD are also present. Significant ternary interaction effects such as ADF, BDF, CEF, and CDE in the ANOVA summary (Table ) indicate that two process variables may also interact with the mixture variables.

Model Coefficients and Interaction Plot

The coefficients of the significant model terms in the total adsorbent loading (qTotal) response are given in eq . This model is expressed in terms of coded variables (with the coding indicated by the apostrophe symbol). The numeric process factors were scaled between −1 and 1 using the respective high (pvhigh) and low levels (pvlow) of each factor as shown in eq .

Concentrations of different mixtures were scaled between 0 and 1 by dividing with the total initial concentration viz. 700 mg/L. The categoric variable F was treated in terms of a two-dimensional vector to denote the three different carbons as recommended by the Design Expert software. The ranges, coding equation, and various levels for each factor are summarized in Table S4 in Supporting Information. The coefficients of the model eq are plotted as Figure to compare across the three carbons and highlight the interaction effects.

Figure 3

Coefficients of the coded qtotal model for the three carbons UN, AT, and MAT.

These coded coefficients (A′–E′) presented in Figure may be compared across the carbons as well as within a carbon. The improvement due to modification is reflected in the increasing positive linear coefficients (A′–C′) in the order of UN < AT < MAT for all three solutes. The coefficients of some binary and ternary interaction terms remain constant for all three carbons. For example, AB, AC, and BC coefficients are constant across carbons, owing to the absence of interaction with the adsorbent. Terms A′B′F′, A′C′F′, and B′C′F′ terms are absent in the model equation (eq ).

Hence, binary interactions between these solutes are unaffected by adsorbent modification. However, other terms highlight the effect of AC modification. For instance, interactions of BTA with pH (B′D′) and ACT with pH (A′D′) depend on the type of carbon.

Also, the coded equation provides insights on the adsorption performance as we increase or decrease a process variable from the design center. At the design center, the coded values of the process variables are all zero. Hence, at neutral pH and mid-dosage (i.e. D′ = 0 and E′ = 0), the first three sets of bar columns of Figure indicate that the single solute preference of all three carbons is as follows: ACT (A) ≈ CAF (B) < BTA (C). However, changing the process conditions disrupts this preference, implying that the single solute preference may be manipulated by altering the process conditions. For example, at pH 3, model predictions revealed that ACT ≈ BTA ≈ CAF, while at pH 10, there is no pattern in the preference order as it depends on the carbon and its dosage. The quadratic effect of pH (D2) is exhibited only when BTA is present, and the terms B′D′ and B′D′2 have negative coefficients, indicating that increasing the pH beyond the center point reduces the BTA loading.

Grouping the terms that have two mixture components, we observe multiple interactions among the solutes as well as with the process conditions. For instance, the binary interaction between ACT and BTA as well as between BTA and CAF is affected by pH, giving rise to differences in the ternary interactions A′B′D′ and B′C′D′. However, since the coefficient of B′C′ is higher than that of B′C′D′, even though pH affects their overall combination, the B′C′ interaction will be positive in magnitude. This indicates that BTA and CAF interact synergistically at any process condition and is likely to exhibit a convex up profile along the BTA–CAF edge of the response surface plot (Figure ). On the other hand, we observe that A′B′ has a lower coefficient value than A′B′D′ (Figure ). Thus, the overall value of A′B′ is positive and negative, when D′ is 1 and −1, respectively. Hence, it is possible that for some pH values lower than the center (i.e. < 6.5), that is, when D′ takes negative values, the overall coefficient can take negative values. Hence, ACT–BTA interactions may be synergetic or antagonistic depending on the pH values. Between pH 3 and 3.95 (D′ = −1 to −0.7), the coefficient of A′B′ is negative, beyond which it is positive indicating antagonistic and synergetic interactions, respectively, for all three carbons. The adsorbent dosage (E′) mainly affects CAF adsorption; while it has a slight influence on BTA, it does not influence ACT adsorption at all. The coefficient values are low, indicating a negligible effect due to dosage.

Figure 4

Response surface of qtotal at various pH values and for different carbon types: UN (a–c), AT (d–f), and MAT (g–i), at a dosage of 0.9 g/L.

Adequacy and Validation of Model Equations

The adequacy of the models for different responses may be quantified in terms of R2, adjusted R2, predicted residual error sum of squares (PRESS), predicted R2, and absence of lack of fit p-value. For the qtotal model, these parameters were 0.9666, 0.9531, 10662.31, 0.9364, and 0.5830, respectively. Low PRESS values and high-predicted R2 values indicate that the models developed may be reliably used for making predictions within the problem domain. High-adjusted R2 and high lack of fit p-value (0.5830) imply that the number of parameters used in developing the final model is neither unduly in surplus nor in deficit.

The models chosen after ANOVA analyses were validated with a new set of twenty experiments as tabulated in the Supporting Information section (Table S5) and shown as a parity plot in Figure S4a. Validation was also carried out at optimal conditions as suggested by Design Expert to maximize qtotal, as tabulated in the Supporting Information Table S6. These results indicate an acceptable prediction capability of within ±10%. The parity between experimental data from the original experimental design and model predictions is shown in Figure S4c. The plot of residuals shown in Figure S4d indicates that the residuals (difference between experiment and model predictions) are normally distributed. The residuals, shown for BTA as an example, in Figure S4e, indicates that the residuals are random with no systematic trends.

Thus, eq can be used for plotting the response surfaces and further optimization. The uncoded version of eq , that is, using actual values of the mixture and process variables for each adsorbent, is used in eq when implementing circular constrained and cyclic optimization methods.

Response Surface Plots of the qtotal Model

The qtotal model (eq ) was used to generate response surfaces (Figure ) for three carbons at three different pH values, that is, nine conditions, at a nominal dosage of 90 mg/100 mL. The two-dimensional contour plots and the optimal qtotal values predicted at these conditions are displayed as Figure S5. These two figures form the basis for the following discussions.

In the response surface plots (Figure ), the three edges of the triangle represent the ACT–BTA, BTA–CAF and CAF–ACT binary systems. In the absence of any interaction between the three components, solute loading profiles for binary and ternary systems will be a linear function of mixture compositions leading to straight line edges and flat planes, respectively. The extent of curvature along the edges in Figure indicates the seriousness of interactions between components. In the present study, at a nominal dosage of 90 mg/100 mL, most of the response surfaces shown in Figure exhibit convex-up curves along the binary mixture planes (or faces), indicating synergetic interactions between the two components. However, inhibitory interactions, that is, those with a concave-up profile, were found in all three carbons at pH 3 for the ACT–BTA binary system alone (Figure a,d,g). This behavior is also evident from the negative coefficients of the A′B′ interaction terms, when the pH value was between 3 and 3.95 (Figure ).

At pH 3, among the three carbons, UN showed the lowest total adsorption capacity. For UN, the synergetic interaction between CAF and BTA led to a maximum qtotal of about 259 mg/g for their binary makeup, while the ACT–CAF binary mixture showed a local maximum of around 238 mg/g. For AT carbon, the performance improved and the maximum loading was around 300 mg/g at the ACT vertex. MAT carbon shows considerably enhanced adsorption (the maximum value for the BTA–CAF binary mixture is about 346 mg/g).

However, this optimization analysis using Figure provides only a single-point prediction at specified process conditions of pH, adsorbent dosage, and adsorbent type. These optima could be local and hence are of limited utility.

Mechanistic Explanation for Adsorption Trends of Different Solutes

Mechanisms such as hydrophobic and π–π interactions are present in the entire pH range,[40] while electrostatic interactions become inhibitive only when the solute species and adsorbent’s surface charges are the same. This phenomenon happens in a narrow pH range determined by the interrelationship between solution pH, pHpzc, and pKa.[41−43] The effect of pH on ACT adsorption can give insight on the adsorption mechanisms of different carbons. In all cases, a low pH was preferred for ACT adsorption. Acid treatment improved ACT adsorption at low pH but not at high pH (Figure ). MAT was found to have noticeably improved ACT adsorption at high solution pH. Wong et al.[44] studied ACT (C0 = 10 mg/L) removal by AC (pHpzc = 2) at a high adsorbent dosage of 2 g/L and reported constant adsorbent loading from pH 3 to 8 and inhibited adsorption at pH 11. They reasoned that as the pKa of ACT was 9.38, it started deprotonating beyond this pH value, and the resulting negatively charged species might have been repelled by the negatively charged AC at higher pH values.

Figure 5

Effect of pH and adsorbent dosage on single-component adsorption.

For a municipal solid waste-based activated biochar (pHpzc = 10) with a dosage of 2 g/L, Sumalinog et al.[43] observed a monotonic decline of adsorbent loading and removal percentage of ACT (C0 = 500 mg/L) throughout the pH range of 2–12, similar to our study.

Hence, when a low solute concentration was used, the dependence manifested only at high pH but at higher concentrations, the dependence was visible throughout the pH range. In Figure , BTA shows a quadratic behavior, indicating the occurrence of maximum adsorption at an optimal pH. The optimal pH values at both adsorbent dosages were 4.8, 5.2, and 3.6 for UN, AT, and MAT, respectively. Here, hydrophobic interactions could dominate BTA adsorption,[42] while electrostatic interactions might hinder BTA adsorption at low and high pH levels,[40] as elaborated below.

This quadratic behavior with varying pH was observed by Sarker et al.[40] for a metal azolate framework (MAF) adsorbent separating BTA from aqueous solutions. The pHpzc for the MAF adsorbent was 8.2. BTA exists in its protonated BTA (+) form at highly acidic conditions (pH < 1.6). This BTA (+) starts to dissociate into neutral BTA and a proton, when the solution pH exceeds its pKa value of 1.6. Thus, at pH 3, about 95% of BTA exist in a neutral form and the remaining as protonated BTA. Between pH 3 and 8.6, neutral BTA dominates. However, neutral BTA dissociates to form deprotonated BTA (−) species, when the solution pH increases beyond its pKa value of 8.6. Thus, the negative species takes over at pH 10. The adsorbents have a positive and a negative surface, when the solution pH is below and above their respective pHpzc, respectively. Here, all three carbons exhibit positive and negative surface charges at pH 3 and 10, respectively. Thus, electrostatic interactions hinder BTA (+) and BTA (−) at pH 3 and 10, respectively, leading to a parabolic profile for BTA adsorption.

CAF exists majorly in its neutral form between pH 3 and 10 since the pKa of protonated and neutral CAF are 0.6 and 10.4, respectively. CAF has a dipole moment due to which dipole–dipole interactions can be present.[45] Moreover, interactions such as H-bonding and π–π interactions may also be relevant. The positive charge on the N atom may interact with a functional group that is negatively charged, and the π electrons present in the 2-nitrophenol ring of CAF may interact with the π electron-rich basal rings of the carbon surface.[46]Figure shows that, at a higher dosage, CAF adsorption is almost pH-independent, whereas at a low adsorbent dosage, CAF adsorption is inhibited with an increase in pH. At higher dosages, the π–π interaction with the basal planes of AC may dominate over other mechanisms. The reduced adsorption at lower dosage and higher pH could be due to the competition between hydroxide ions and CAF for rarer adsorption sites (Ravi et al. 2020), even though the N atom of CAF may interact with the negatively charged carbon surface.

For AC fibers having a pHpzc of 2.8, Beltrame et al.[41] have reported reduced CAF adsorption when pH was beyond 7. They had used 1 g/L of AC fiber for the removal of 500 mg/L CAF and attributed the decrease to electrostatic repulsions. In contrast, Portinho et al.[47] observed no dependence of pH on CAF adsorption when 1 g/L of grape stalk AC was utilized to remove 20 mg/L of CAF. However, this pH independence could be due to surplus of adsorbent exposed to relatively low initial concentrations, as was the case with ACT,[44] which was discussed above. The above two studies support the present observation, that is, at a high dosage or low initial concentration, the pH dependence of CAF adsorption is almost non-existent.

The discussions given above demonstrate the ability of MPV design model predictions to aid holistic investigation of adsorption performance. It also enables additional insights on the adsorption process, and the conclusions that are made find support from earlier literature studies.

Optimization of Model qtotal

The model developed in Section (eq ) can be analyzed and used for determining the global maximum in qtotal. Further, the evolution of local optima to the maximum qtotal along the path of steepest ascent may be described in the ternary composition space using ridge analysis.[31] When tracking these local optimal solutions, the process variables may or may not be kept fixed as explained below.

Conventional Ridge Analysis for Mixture Design

First, the settings of both the mixture and process variables that led to highest adsorbent loading is determined using eq for each adsorbent. This is the reference condition with respect to which other local optimal solutions are compared. The global maxima conditions and qtotal values of eq in the non-coded form were obtained using particle swarm optimization routine of MATLAB and are tabulated in Table .

Table 3

Conditions That Result in Globally Optimal qtotal Value, Considering pH as the Continuous Numerical Variable between 3 and 10

	C0,ACT mg/L	C0,BTA mg/L	C0,CAF mg/L	pH	dosage mg/100 mL	qtotal mg/g
UN	0	530.77	169.23	4.91	60	279.69
AT	527.15	0	172.85	3	120	303.72
MAT	0	387.73	312.27	4.09	60	368.53

Now, to implement the conventional ridge analysis, the process variables must be fixed. Fixing the pH and dosage at values given in Table converted the higher-order MPV model (eq 2) to a quadratic mixture design model (eq ). This enabled the application of ridge analysis theory detailed in Section .

The pH and dosage corresponding to global optimum (Table ) were specified for each carbon as they would be the most logical choices for a priori specification. In Figure , the contour plots and the ridge analysis results, plotted as hollow symbols, are shown for each carbon. From the contour plots, it can be seen that MAT AC is distinctly superior to the other carbons in terms of total adsorption. The contour values show higher values of qtotal for MAT carbon. The topology of AT carbon is drastically different when compared to those of UN and MAT carbons.

Figure 6

Loci predicted by conventional ridge analysis (a–c) and circular constraint optimization (d–f) at pre-specified globally optimal process variables. Locus of solutions from the conventional ridge analyses (a–c) terminates once it reaches the binary edge of the triangle. It reaches global optimum only for a specially chosen point P (a–c). The circular constraint method reaches the global optimum by moving along the binary edge. (P: specially chosen point, A: an arbitrarily chosen point, C: centroid, and G: global optimum co-ordinate).

When the centroid C (700/3, 700/3, and 700/3) in the triangular domain was chosen as the focus and the conventional ridge analysis method outlined in Section was used, the locus of optimal compositions represented by hollow circles did not pass through or terminate at the global optimum, G. However, to ensure that the ridge analyses path shown by hollow triangles could be made to pass through the global optimum, a suitable focal coordinate, P, could be identified for all the three carbons. This coordinate P was located at the point of intersection of the binary edge and the perpendicular drawn from the global optimum coordinates G. The response spaces traversed by circles drawn from this coordinate as a center could get into the higher response region and eventually even the conventional ridge analysis (hollow triangles in Figure a–c) could terminate at the global optimum G. With focus at P, it is obvious that one of the circles will intersect the triangle’s other edge tangentially at the global optimum mixture composition. Since the process variables were fixed a priori at global optimum conditions, the mixture variables should also converge to the global optimum settings identified in Table . There are three limitations to this conventional ridge analysis method though it enjoys the advantage of theory as described below. Theory of conventional ridge analysis is given in the Supporting Information section.

The steepest ascent trajectory of maximum total adsorption may not always terminate or pass through the global optimum because of the inappropriate selection of the starting point, that is, the focus (Figure a–c).

Ridge analysis locates the maximum total loading on concentric circles (Figure a). When the circle projects out of the mixture design space, infeasible optimal solutions such as negative compositions may arise on those circles.

The pH and dosage process variables set a priori at global optimum (Table ) may actually turn out to be sub-optimal for mixture compositions along the ridge analysis trajectory except at the globally optimal ones.

This analysis, despite its limitations listed, reveals the preference of AT carbon toward ACT and that both UN and MAT show similar trends of preferring the BTA–CAF binary mixture. Furthermore, this analysis would serve to indicate optimal feed concentrations that maximize adsorbent utilization for specified process variables, be it the globally optimal process conditions or those dictated by operating plant constraints.

Improved Optimization of qtotal with Process Variables Fixed at Their Global Optima

In this analysis, the process variables such as pH and dosage are still fixed at the globally optimal values. The first and second limitations (a,b) of ridge analysis as listed above in Section were overcome by using the composition constraints (eq ) in addition to the mixture constraint (eq ). These constraints were implemented using MATLAB routine fmincon, which was first validated as reported in Table S7 of Supporting Information. The constraint to lie within the composition domain is schematically illustrated in Figure b. After addition of the composition constraints (eq ), it was observed that the steepest ascent trajectory could be made to go through the global optimal values irrespective of the focal values, A or C, chosen as starting points (solid line and dots in Figure d–f).

When centroid C is the focus, comparing both ridge analysis (hollow circles) and constrained optimization (solid line), we observe that the conventional ridge analysis reached the triangle’s edge but it failed to progress toward the actual global composition present on the binary edge. This was because of limitation “b” in the list mentioned in the previous Section . However, the locus of compositions predicted by constrained optimization not only reached the binary edge but could move along it and reach the global optimum in both cases such as with focus as centroid and an arbitrarily chosen point, A.

However, the ridge analysis so far had pre-specified process conditions corresponding to the global optimum. However, as per limitation c in Section , these process conditions may actually turn out to be sub-optimal at other compositions. Piepel et al.,[48] observed that a two-step procedure involving optimization of the mixture composition at specified process conditions and then optimizing the response surface of the process conditions may lead to misidentifying the design space as well as sub-optimal settings of the mixture and process variables. They attest that when interactions arise between mixture and process variables, their optimal settings will be mutually dependent on each other. To identify optimum values of process and mixture variables simultaneously, when maximizing qtotal, the ridge analysis for mixture design was extended to the MPV design.

Optimization of qtotal Considering Mixture and Process Variables

In the cyclic optimization method outlined in Section , both process and mixture variables are allowed to change simultaneously. This facilitates the comprehensive and simultaneous investigation of the entire MPV design space. An additional advantage is that the maximum adsorbent loadings along the evolving path may potentially be higher than those obtained by the previous two strategies (Sections and 3.2.2) because pH and dosage in addition to composition are getting adjusted during the search for the optimum. In this analysis, the choices of foci were based on two criteria. First, the focal points should be far from the global optimum for the chosen carbon, so that the steepest ascent path can be sufficiently long and distinct. The second criterion is that the focus may be either the centroid or an arbitrarily chosen point subject to criterion 1. The results from these analyses are plotted as Figure for an arbitrary point and as Supporting Information Figure S7 for the centroid. The optimum responses for qtotal along the path of steepest ascent are shown in Figure a. These optimal responses are labeled as “MPV”. Also shown as comparison are the trends where only mixture compositions were optimized after setting process variables at their global optimum values. These are labeled as “Pure Mix”. In Figure b–d, plots for loci of pH, individual component concentrations, and the steepest ascent path in the ternary composition space corresponding to the locus of maxima in qtotal are presented. In Figure b, the global optimum in pH is also shown for comparison. The optimum dosages in MPV optimization, however, did not change from the global optimum values, reported in Table , for the three carbons and hence the plots for dosage are not shown.

Figure 7

Evolution of (a) locally optimal responses and their corresponding (b) process and (c) mixture variables with increasing radial distance from the focus. (d) Loci of mixture compositions along the path of steepest ascent during cyclic optimization with the MPV model and its composition and circular constraints.

For UN, at initial stages, pH 3 was found more suitable than the globally identified pH value of 4.91 as a higher qtotal can be achieved. This is shown in Figure b. In the case of MAT (Figure a), the locus of maximum qtotal from mixture composition optimization (“Pure Mix”) nearly overlaps with that of MPV optimization. The optimal pH over the entire path varies between 4 and 4.5 (Figure b), and this range represents a small deviation from the global optimum pH of 4.09. For the UN adsorbent, ACT and CAF had similar optimal composition trends (Figure c). For AT, the ACT and BTA trends were opposite while that of CAF was ambivalent within a relatively narrower composition range. For MAT, however, BTA feed compositions were relatively invariant when compared to those of ACT and CAF with the latter two exhibiting opposing trends. The optimal pH trend seemed to correlate more closely with the trend of BTA in both UN and MAT. For AT, the optimal pH trend was more closely correlated with CAF. These plots (Figure b,c) are also useful in defining feed concentration ranges where the adsorption process is either independent of or most sensitive to pH.

Due to the circular domain constraint, the variables are forced to lie on the circle rather than inside the circle. Hence, the evolving response as well as process variables need not be constant or change monotonically with radius as there is a higher or lower local optimum relative to the previously identified one along the trajectory originating from the focus. These non-monotonic optima are revealed as humps in the plot of qtotal versus radial distance (Figure a). However, among the carbons considered, only AT is the exception as it displayed a non-monotonic behavior toward the global optima. As depicted in Figure a, when the focus is at the BTA–CAF edge, a hump indicating local optima is observed at an intermediate location for AT carbon (red lines) with pH nearing 6.5. Here, qtotal is about 298 mg/g. With a further increase in the radial distance, a slightly better optimum (qtotal of 303.72 mg/g) involving ACT and CAF is identified at pH 3. None of the composition triangles has contour plots embedded in them because the process conditions keep changing along the trajectory originating from the focus.

MPV Optimization for Percentage Removal

The removal of solutes from the aqueous stream is equally an important objective as maximizing the utilization of typically expensive adsorbents.[49−51] Furthermore, the two objectives are not equivalent. Hence, percentage removal was also modeled from the experimental data using Design Expert 11 (Stat-Ease, Inc., Minneapolis) and subsequently optimized. The total percentage removal (PRtotal) is defined as follows:where C0,total is the total initial concentration and Ce,total is the final equilibrium concentration in the solution.

The mixture and process variables optimized for maximum adsorbent loading may not be optimal for maximum total percentage removal as well. By definition, the percentage removal is based on mass of the solute adsorbed per unit mass of solute in the feed, while the adsorbent loading is based on the mass of solute adsorbed per unit mass of adsorbent. The loading of the solute per unit mass of the adsorbent may be maximized by either increasing the amount adsorbed or minimizing the dosage. On the other hand, the percentage removal may be maximized by increasing the adsorbent dosage and/or lowering the equilibrium concentration of the solute (Ce) by using a better adsorbent.

The coded form of PRtotal given as eq was found from experimental data through regression analysis and contains only the statistically significant terms. Here, mixture compositions vary from 0 to 1 and process factors vary from −1 to +1 as reported in Supporting Information Table S4. A parity plot of the predicted and actual % removals is also given in Supporting Information Figure S4b, which indicates an adequate fit by the model. Here, validation data set that was independent of the original design was used (Table S5). This model (eq ) had high R2 (0.982), adjusted R2 (0.976), and insignificant lack of fit (p-value = 0.2248). The ANOVA for PRtotal is given as Table S9 in Supporting Information.

The cyclic optimization technique with centroid as the focus, involving both mixture and process variables was carried out to maximize PRtotal. The results are depicted in Figure . The locally optimal total percentage removal values along the steepest ascent path increased with an increase in the radial distance from the focus (Figure a). Also compared in Figure a in dashed lines is the PRtotal trend when qtotal was maximized instead.

Figure 8

Evolution of PRtotal (a), qtotal (b), pH (c), and mixture variables (d) when PRtotal and qtotal were individually maximized using the cyclic optimization method. The solid line indicates MPV maximization of PRtotal, while the dashed line indicates MPV maximization of qtotal. The focus is at the centroid.

The percentage removals here were found to be much lower. Similarly, the qtotal values were lower when PRtotal is maximized. This is shown in Figure b. For example, when PRtotal was maximized for UN, the locally optimal PRtotal values were between 40 and 46% (Figure a, solid line) and the local optimum values of qtotal values were in the range 247–250 mg/g with a maximum in between (Figure b, solid line). When qtotal was maximized, the local optimum values of PRtotal values fell in the range 20–25% (Figure a, dashed line). On the other hand, the local optimum values of qtotal increased substantially to the range 254–280 mg/g (Figure b, dashed lines). Similar observations can be made for other carbons.

When the total adsorption capacity was maximized, qtotal increased with the radius but the corresponding PRtotal monotonically decreased for AT carbon but decreased only after a while for UN and MAT carbons (Figure b). Thus, for a certain combination of mixture composition and process variables, a compromise solution that achieves the highest possible values of both qtotal and PRtotal will be useful for the water treatment plant. Even when the composition had all the solutes, (near the centroid) that is, R ≈ 0, the optimal pH was between 4 and 5 for all the carbons and gradually reached the respective globally optimal pH (Figure c). In Figure d, we observe how the feed compositions initiated from the centroid for different adsorbent change along the respective steepest ascent paths.

These steepest ascent paths are different for individual maximization of PRtotal and qtotal, again illustrating the stark difference between the two optimization strategies. These paths represent the loci of the local optima before eventually terminating at the respective global optima. The global optima for PRtotal are reported in Table . As indicated in Table , AT and MAT have higher BET surface areas and hence adsorb better than UN.

Table 4

Conditions That Result in Globally Optimal PRTotal Value, Considering pH as a Continuous Numerical Variable between 3 and 10

AC	C0,ACT mg/L	C0,BTA mg/L	C0,CAF mg/L	pH	dose mg/100 mL	% removal
UN	0	700	0	5.0	120	45.60
AT	0	700	0	4.1	120	58.75
MAT	0	700	0	4.7	120	58.00

Irrespective of the adsorbent, the steepest ascent path based on PRtotal terminated at the global optimum located at the vertex of the ternary diagram corresponding to a pure BTA. Hence BTA alone is the preferred feed for maximum PRtotal under acidic conditions. This indicates that when there is sufficient amount of adsorbent for all the solutes to adsorb, BTA is preferred by all the carbons. The binary mixtures are preferred at low dosage values, when there is competition for active sites. The dosage corresponding to the maximum percentage removal hits the upper bound for all three carbons and hence is not reported in Figure . This is expected because, generally, a higher dosage leads to a higher removal.

It is to be noted that several such versions of Figure , depicting the locus of optimal conditions and path of steepest ascent may be easily generated when the feed composition originates from a focus from anywhere else in the triangle. It is also observed in Figure a,b that when either qtotal or PRtotal is maximized individually, the other response is clearly adversely affected. Hence, both responses may be optimized together to arrive at the best possible compromise solution.

MPV Optimization Combining Both qtotal and PRtotal

In an operating plant, it is preferable to maximize the utilization of the adsorbent as well as ensuring high removal of solutes from the aqueous feed solution. In other words, both qTotal and PRtotal have to be maximized and one way of accomplishing this is to combine the model equations developed for both criteria (eqs and 9) in a suitable manner into a single objective function. Due to difference in units, the combined objective function cannot be directly expressed in terms of a sum of qtotal and PRtotal. Furthermore, the numerical values of these responses are of different magnitudes, 0–59% for PRtotal and 0–368 mg/g for qtotal. Freitas[52] recommended that if the magnitudes of the responses varied widely, it was preferable to scale them. The qtotal and PRtotal responses were scaled by their global optimum values, respectively. This scaling makes both total adsorbent loading and total percentage removal dimensionless and forces them to range from 0 to 1. The two scaled responses were combined into a single weighted objective function as shown by eq .

The weights are indicative of the relative importance of the responses. The weights are usually normalized so that they add up to unity.[53] In the present case, different combinations of weights attached to qtotal and PRtotal were tried for the three carbons as shown in Figure . PRtotal decreased rapidly when the weight for qtotal was increased beyond a threshold value. This knee junction forms the best compromise of having high qtotal and PRtotal values concurrently. The individual maxima and the best compromise weight values are summarized in Table S10 for the three carbons.

Figure 9

Fronts of feasible solutions ranging from global maximum for PRtotal and global maximum for qtotal.

The weights (w1 and w2) that led to highest qtotal and highest PRtotal were obviously (1, 0) and (0, 1), irrespective of the carbon. Even though many other desirable outcomes may be stipulated, the most common one is the combination of weights that finally leads to high value of both the responses. Table S10 indicates that higher weight has to be attached to qtotal even after scaling the response equations with qtotal,max and PRtotal,max. Equal weight combination (i.e., 0.5 and 0.5) did not lead to high values of both qtotal and PRtotal, as evident in Figure . The ε-constraint optimization method where one of the objective functions is provided as a constraint was also tried but the compromise solutions were inferior to those obtained from the weighted sum optimization method.

Utility and Importance of the Proposed MPV DOE and Optimization Approach

The MPV design provides an experimentally economical method to study the multicomponent adsorption performance involving both mixture and process factors. Industrial adsorbents are often expensive, and the method developed in this work would eliminate costly and inaccurate trial and error in adsorbent selection and identification of optimal conditions for the encountered combination of pollutants.

The significant advantage of using the novel cyclic optimization technique is that from any starting point in the feasible composition space, one may trace the steepest evolution of the local optima en route to the global optimum conditions. This additional knowledge may be valuable in a practical scenario where the global optimal condition may not be viable. For instance, in a wastewater treatment plant dealing with three pollutants, global optima involving a binary mixture or a single component would not be of any use but one of the local optima identified that involves all the three components in the feed could be actually implemented. In the present work, it appears that none of the three carbons prefers adsorption of the ternary mixture. Thus, if a ternary mixture has to be treated in the wastewater plant, the process has to be operated at one of the locally optimal conditions only. For a given ternary feed mixture, a favorable local optimum involving a moderate pH, a suitable AC, and a low adsorbent dosage may be identified and would be of practical utility.

In a preliminary study, Retnam[54] has confirmed that even for continuous multicomponent adsorption in a packed column, MAT carbon performed better than UN and AT carbons at the optimal feed compositions and pH values suggested by the qtotal model used in this work. However, more analyses would be required to extend this work to continuous systems. Another advantage in the developed method is that both process variables and mixture variables are accounted for in the adsorbent selection.

Conclusions

The MPV design coupled with constrained optimization facilitates economical and holistic investigation of multicomponent adsorption. The distance-based optimality criterion dispersed the selected points uniformly throughout the design space (Supporting Information Figure S1). Statistical significance of the main factors and interactions involving the mixture and process variables (pH, adsorbent dosage, and type of adsorbent) was captured using ANOVA (Table ). The statistical analyses revealed that many interactions between factors could be more influential and significant (p-values <0.001) than the main factors alone (p-value <0.01). The model suggests a strong interaction between the process and mixture variables for example, pH and BTA concentration (p-value < 0.0001). Furthermore, the mixture variables interacted significantly among themselves as well, for example, CAF with ACT and BTA (p-values < 0.0001). The total adsorbent loading and percentage removal responses were succinctly described in terms of robust empirical models (eqs and 9). These models for qtotal and PRtotal have good predictive capability for single, binary, and ternary systems (R2 > 0.96, Adj. R2 > 0.95, model prediction R2 > 0.93) and insignificant lack of fit (p-value > 0.22).

The present study demonstrates the potential of microwave and acid treatments in considerably enhancing the performance of adsorbents. The qtotal of AT and MAT adsorbents improved by 8.6 and 31.8% relative to the unmodified (UN) adsorbent, respectively. The total percentage removal of both AT and MAT adsorbents improved by about 28% over the UN adsorbent.

The application of ridge analysis may lead to the identification of infeasible solutions lying outside the allowed composition space. Furthermore, the path of steepest ascent describing the evolution of locally optimal responses (qtotal and PRtotal) did not always eventually culminate in the global optimum (Figure a,b). Even better solutions may have been obtained along the steepest ascent path had the process variables also been varied in addition to mixture variables. Invoking the circular and individual composition constraints enabled the attainment of feasible solutions within the mixture-composition space (Figure d–f). The circular constraint optimization method could involve both mixture and process variables, but the steepest ascent path did not reach the global optimum. This limitation was overcome by the novel cyclic optimization technique (Figure ).

The adsorbent utilization and percentage removal were optimized individually as well as together. When qtotal and PRtotal were maximized individually, the corresponding PRtotal and qtotal however decreased (Figure ). Hence, it is better to optimize both qtotal and PRtotal together for better utilization of the adsorbent, which are usually expensive while at the same time not compromising with purification of the wastewater (Figure ). The extension of the enhanced ridge analysis-based optimization methods developed in this work for batch adsorption systems to continuous multicomponent systems is recommended. It is also recommended to understand the dynamic competition between different solutes during adsorption through multicomponent batch kinetic studies.