Predicting adsorption of water/organic mixtures using molecular simulation

Predicting Adsorption of Water/Organic Mixtures Using Molecular Simulation Miguel Jorge and Nigel A. Seaton School of Engineering and Electronics, University of Edinburgh, Edinburgh, EH9 3JL, U.K. The use of Monte Carlo simulation to predict the adsorption of mixtures of polar and nonpolar species on acti®ated carbon was in®estigated using water and ethane on BPL carbon as a prototype system. The structure of the adsorbent was modeled by an array of slit-shaped pores, characterized by a pore-size distribution. The chemical heterogeneity of the carbon was taken into account by including oxygen-containing sites on the surface of the pores. The pore-size distribution was obtained from pure-ethane adsorption on the same carbon sample, while the concentration and distribution of surface sites were determined by analyzing pure-water adsorption. Model predictions agree well with experimental multicomponent data. Introduction A substantial proportion of industrial adsorption processes Ž involve carbon adsorbents Bansal et al., 1988; Sircar et al., .Ž 1996 . In many important applications such as air purifica- . tion and personal protection , there is a need to selectively adsorb organic species from gas streams containing a certain amount of water vapor. It is well known that the presence of water in the feedstream can significantly affect the efficiency of the separation process. Reduced adsorption capacities and earlier breakthrough times have been reported for adsorption of a wide range of organics on activated carbon in the Ž presence of water Gong and Keener, 1993; Scamehorn, 1979; . Shin et al., 2002 . The development of a method to predict the adsorption of mixtures of water and organic species on activated carbon is, therefore, essential for the efficient design of such processes. This in turn requires a detailed understanding of the fundamental aspects that affect adsorption of both polar and nonpolar components. The most common approach for predicting multicomponent adsorption is based on classical thermodynamic methŽ. ods. Of these, the ideal adsorbed solution theory IAST Ž. Myers and Prausnitz, 1965 is the standard method used in industry. It is based on the concept of spreading pressure and assumes that the adsorbed phase is an ideal solution in equilibrium with the bulk gas phase. For systems such as mixtures of water and hydrocarbons adsorbed on activated carbon, the IAST is likely to perform badly, since this system is far from ideal. Other thermodynamic methods have therefore been Correspondence concerning this article should be addressed to N. A. Seaton. Ž. proposed. To our knowledge, Okazaki et al. 1978 were the first to attempt a prediction of multicomponent adsorption in this type of system. In their method, water is adsorbed by capillary condensation, while the adsorption of the organic species is the result of three contributions: adsorption on the dry adsorbate surface, dissolution into the condensed phase, and adsorption onto the wet surface. They applied this method to binary adsorption of water and several organic Ž. species acetone, methanol, benzene, and toluene on activated carbon. Although the predictions were in reasonable agreement with experimental data, the method has the significant disadvantage of requiring large amounts of data as inŽ put apart from the pure-component adsorption isotherms, it also requires liquid-phase isotherms for the organic species . in water and vapor᎐liquid equilibrium data . A later method, Ž. based on potential theory, is due to Doong and Yang 1987 . Their method uses the concept of maximum available pore volume, together with the Dubinin-Radushkevich equation Žand later, the Dubinin-Astakhov equation; Doong and Yang, . 1988 with parameters obtained from the pure-component isotherms. A system with no lateral interactions between adsorbed species is assumed, but the maximum available micropore volume for a given species is reduced by the amount adsorbed of all other components. This method was seen to perform well for the same systems studied by Okazaki et al. Ž. 1978 . A similar approach was used by Lavanchy and Stoeckli Ž. 1999 to predict binary adsorption of water and 2-chloropropane on activated carbon. Their method is essentially a simplified form of the Doong-Yang method using the DubininAugust 2003 Vol. 49, No. 8AIChE Journal 2059 Astakhov equation. Lavanchy and Stoeckli also applied IAST, coupled with the Dubinin-Astakhov equation to represent the pure-component data, to adsorption of waterrorganic mixtures. They concluded that IAST performs reasonably well when the organic component is soluble in water, while their implementation of the Doong-Yang method is required for water-insoluble organics. Similar conclusions were drawn Ž. more recently by Linders et al. 2001 . Finally, the efforts of Ž. Appel et al. 1998 to describe multicomponent adsorption by using mixture virial coefficients are noteworthy. These authors successfully applied this theory to binary waterrhexane adsorption on activated carbon, by using the equation of Talu Ž. and Meunier 1996 for the pure-component water isotherms. However, the calculation of the virial mixture coefficients requires the fitting of multicomponent data, and, hence, their method was only applied in a correlative manner. The thermodynamic methods described previously have some important disadvantages: they become inaccurate as the Ž systems move away from ideality in the sense of Raoult’s . law , due to chemical dissimilarity among the adsorptive species; they require a significant amount of experimental data; and they give us no information on the structure of the adsorbent. Recently, an alternative approach, based on statistical mechanical methods such as density functional theory and molecular simulation, has begun to be applied to multiŽ component adsorption prediction see Razmus and Hall, 1991; Segarra and Glandt, 1994; Gusev and O’Brien, 1998; Pan et al., 1999; Vuong and Monson, 1999; Macedonia and . Maginn, 1999; and references therein . In particular, molecular simulation methods have been shown to be extremely powerful in predicting multicomponent adsorption in activated carbon as well as in characterizing the adsorbent itself Ž. Davies and Seaton, 2000 . Molecular simulation can be used to obtain information about the chemical and physical characteristics of the adsorbent that affect the adsorption of the molecular species of interest, using a limited amount of experimental data as input. With this information, a model adsorbent, representative of the real adsorbent, can be constructed and used to predict adsorption in the real adsorbent. The effect of physical and chemical changes on adsorption can also be investigated using the model adsorbent, thereby using molecular simulation as an adsorbent design tool. In this article, we attempt to develop such a model for the case of waterrorganic mixtures. Models for Activated Carbon Activated carbons are composed of graphite plates stacked Ž. together in a more or less irregular way McEnaney, 1988 . The interstices of these plates form a network of interconnected pores, between which adsorption occurs. These carbons are usually microporous and have a very strong affinity toward hydrocarbons. However, in spite of their generally hydrophobic nature, most activated carbons used in industrial applications are seen to strongly adsorb water. This effect is due to the presence of polar groups that appear on the edges of the graphitic plates or on defect sites, as a consequence of the activation process. The most common types are oxygencontaining groups, such as carboxyl, carbonyl, and hydroxyl Ž. Boehm, 1994 . These groups form nucleation sites for water, increasing the hydrophilic properties of the carbon. Thus, the structure of the carbon has a strong effect on adsorption of all chemical species, while the chemically heterogeneous sites affect mainly the adsorption of polar species. In the case of waterrorganic mixtures, both types of heterogeneity must be taken into account when developing a model for the adsorbent. The most widely used model of the structure of activated carbon is an array of slit-shaped pores, characterized by a Ž. pore-size distribution PSD . This model takes advantage of the shape of the graphite platelets and assumes that the entire carbon can be represented by a bundle of infinite pores of slit shape, which are in equilibrium with the bulk gas phase. The total adsorption on the carbon can then be related to the Ž.Ž PSD via the adsorption integral equation AIE Seaton et . al., 1989 ⬁ NT,P,ysfw ␳ T,P,y,wdw 1 Ž. Ž. Ž. Ž . H ii 0 Ž. where NT,P,yis the total adsorption of component iat i temperature T, bulk-phase pressure P, and composition y; i Ž. ␳ T,P,y,wis the molar density of the adsorbent at the i Ž.Ž . same conditions in a pore of width w;fw sdV rdw is mic the PSD; and Vis the specific micropore volume of the mic adsorbent. The single-pore isotherms can be calculated in a number of different ways, ranging from macroscopic models like the Kelvin equation to the more complex molecular simulation methods. Attempts to include chemical heterogeneity in model carŽ. bons have so far been few. Segarra and Glandt 1994 have modeled activated carbons as a disordered stacking of diskshaped graphitic platelets. Surface polarity was represented by dipole moments smeared along the edges of the plates. Simulated adsorption isotherms in this model were compared to experimental data, but the agreement was only qualitative. Ž. Subsequent work Gordon and Glandt, 1997 suggests that the smeared dipoles are too homogeneous to represent the polarity of real adsorbents, and that a discrete distribution of polarity would be more realistic. A later model, by McCallum Ž. et al. 1999 , combined a discrete representation of the surface sites with a PSD model for the carbon structure. The PSD was determined by using density functional theory to analyze a nitrogen adsorption isotherm, and surface sites were randomly distributed on the surface of the pores. Simplified molecular models, based on square-well potentials to mimic hydrogen bonding, were employed to represent water molecules and surface oxygen atoms. The total concentration of polar sites on the carbon was determined experimentally by titration methods and the interaction potential between water molecules and surface sites was obtained from a fit to low-pressure water adsorption data. Agreement with experiment was reasonable. However, the potential models ignore the long-range nature of the electrostatic interactions among water molecules, and between the water molecules and polar surface groups. This makes it difficult to extend the model to Ž more general cases such as other temperatures, multicompo- . nent adsorption . A more recent model of activated carbon to include chemical heterogeneity was proposed by Brennan Ž. et al. 2001 . The structure of the adsorbent is represented by a stacking of graphitic plates of several sizes, which is optiAugust 2003 Vol. 49, No. 8 AIChE Journal2060 mized to reproduce experimentally determined carbon᎐ carbon radial distribution functions by a procedure termed Ž reverse Monte Carlo simulation Thomson and Gubbins, . 2000 . Chemical heterogeneity was included by randomly placing polar carbonyl sites on the edges of the plates. Although this model provides us with a more realistic representation of activated carbon, it is also extremely complex, which carries disadvantages such as very computationally demanding simulations and difficulty of application to a wide range of adsorbents. To our knowledge, no comparison of this model with an experiment has yet been presented. The model proposed here is based on a PSD to represent the physical structure of the activated carbon, combined with polar surface sites to describe the chemical heterogeneity. This distribution is obtained from the AIE using experimental data for pure-ethane adsorption together with grand Ž. canonical Monte Carlo GCMC simulation to calculate the single-pore isotherms. We have included chemical heterogeneity in the model by placing polar oxygenated sites on the surface of the pores, the concentration and distribution of which are obtained from an analysis of experimental purewater adsorption. This model is expected to capture all the fundamental characteristics of the adsorbent that affect the adsorption of both polar and nonpolar species, while retaining a degree of simplicity that will enable its use in an industrial context. Molecular Simulation in Single Pores We have used GCMC simulation to calculate equilibrium adsorption isotherms in model pores. The Monte Carlo method is a stochastic process based on the laws of statistical mechanics. In effect, one samples a large number of molecular configurations that are consistent with a given set of Ž.Ž. macroscopic thermodynamic variables Hill, 1960 . For most adsorption calculations, the grand canonical ensemble is the most convenient. In this ensemble the temperature, volume, and chemical potential are kept constant, while the total number of molecules is allowed to fluctuate. This situation is analogous to an adsorption experiment, and, thus, adsorption isotherms can be obtained from GCMC simulations by calculating the equilibrium adsorbed density at each set of conditions. In practice, instead of fixing the chemical potential, it is more convenient to specify properties of the bulk phase, Ž. Ž. such as pressure Pand composition x. One must then i employ an appropriate equation of state to calculate the chemical potential of the bulk phase. In the simulations reported in this article, the Peng-Robinson equation of state Ž. Sandler, 1989 was used for this purpose. The GCMC simulation method has been widely used and is well documented ŽŽ. see, for example, Allen and Tildesley 1989 and Frenkel and Ž.. Smit 1996 . Details of the implementation for this particuŽ lar system have been published elsewhere Jorge and Seaton, . 2002a . The simulation cell is rectangular and is bounded in the z-dimension by the pore walls, the distance between which Ž. defines the pore width w. In both directions parallel to the Ž. walls xand y, periodic boundary conditions were used to represent a semi-infinite pore. The length of the cell in these directions was 3 nm, which is sufficient to avoid any finite-size Ž. effects Jorge and Seaton, 2002b . Ethane molecules were Table 1. Geometric and Potential Parameters U for the Molecular Models Used in This Work Parameters Graphite Ethane Water Carbonyl Ž. ␴ nm 0.34 0.3512 0.3166 0.296 y1 Ž. ⑀ Jⴢmol 232.8 1162.3 650.2 879.6 Ž. ⌬nm 0.335 ᎏᎏᎏ y3 Ž. ␳ nm 114 ᎏᎏᎏ sŽ. Bond length nm ᎏ0.2353 0.1 0.1233 Ž. Bond angle ⬚ᎏᎏ109.47 ᎏ yŽ. qe ᎏᎏy0.8476 y0.5 qŽ. qe ᎏᎏq0.4238 q0.5 U Ž. Graphite parameters were taken from Steele 1974 ; ethane parameters Ž. from Cracknell and Nicholson 1994 ; water parameters from Berendsen Ž. et al. 1987 ; carbonyl parameters from Jorgensen and Tirado-Rives Ž. 1988 . Ž. represented by two Lennard-Jones L-J sites, one for each CH group, with parameters taken from the work of Crack3 Ž. nell and Nicholson 1994 . The point-charge SPCrE model of Ž. Berendsen et al. 1987 was chosen to represent water molecules. The interaction between an L-J center in an adsorbate molecule and each of the pore walls is given by the Ž. 10-4-3 potential of Steele 1974 for graphite. The oxygenated sites on the surface of the carbon were represented by carbonyl groups, described by the OPLS potential for C⫽O Ž. groups in aminoacids Jorgensen and Tirado-Rives, 1988 . This model is composed of an L-J center and a negative point charge for the oxygen atom and a positive point charge on the carbon atom in the basal plane of graphite, as described Ž. in a previous publication Jorge et al., 2002 . Geometric and potential parameters for water, ethane, graphite, and carbonyl sites are given in Table 1. The potential between two point charges is long-ranged and its calculation requires the use of special techniques. The full extent of these long-range interactions was accounted for using the method of Heyes Ž. and van Swol 1981 . Cross-species Lennard-Jones parameters were calculated using the Lorentz-Berthelot combining rules. When comparing simulation results with experiments, the former needs to be corrected to account for the difference between absolute and excess adsorption. The result of a simulation is the total adsorption of each species in a model poreᎏabsolute adsorption. On the other hand, adsorption experiments measure the difference between absolute adsorption and the amount of adsorbate that would be present under the same thermodynamic conditions if there were no interactions with the adsorbentᎏexcess adsorption. The simulated absolute adsorption isotherms are converted to excess isotherms using the method proposed by Davies and Seaton Ž. 1999 . Pore-Size Distribution Ž. The PSD is obtained from the AIE Eq. 1 , using experimental data for the adsorption of a pure component, which gives the number of moles adsorbed, as a function of temperŽ. ature and pressure, NT,P,ys1 . As mentioned previously, i the single-pore isotherms were calculated by GCMC simulaŽ. tions. Equation 1 is then solved for the PSD, fw. This is an ill-posed problem, as in principle an infinite number of PSDs can satisfy the experimental data. We have used the method August 2003 Vol. 49, No. 8AIChE Journal 2061 Ž. of Davies et al. 1999 to solve the AIE. A brief outline of this method is given next. First, the AIE is discretized, using a number of quadrature points that is less than or equal to the number of experimental data points. In order to avoid excessive spikiness of the PSD, a regularization factor is incorporated in the minimization procedure. The optimal regularization factor is determined by analyzing the error of the fit, using both generalized cross-validation and L-curves. This procedure eliminates the need to restrict the PSD to a specific functional form, while still ensuring that a physically meaningful result is obtained. When establishing the quadrature intervals, one must take into account the concept of ‘‘window of reliability,’’ first Ž. introduced by Gusev et al. 1997 . Essentially, this establishes lower and upper limits to the pore sizes that can be reliably identified in a PSD analysis. The lower limit is the smallest pore that the adsorbate of interest can enter, while the upper limit is the largest pore for which the single-pore isotherms are still linearly independent. For large pores, the adsorption will occur virtually independently on each pore wall, and the pores will then become indistinguishable from each other. In this method, all the pores above the upper limit of reliability are accounted for in a single quadrature interval. The limits of reliability depend on the adsorbate, on the temperature, Ž. and on the pressure. Davies and Seaton 2000 have recently shown that the window of reliability is extremely important when predicting adsorption, while the resolution of the PSD plays only a minor part. According to the preceding procedure, the discretized AIE becomes m NT,P,yffw ␳ T,P,y,w ␦ w2 Ž. Ž. Ž. Ž . Ý ininn ns1 where mis the number of quadrature intervals used in the analysis, with m-1 intervals lying within the window of reliability and one interval accounting for the adsorption in all the pores above the upper limit. Figure 1. Pore-size distribution of BPL activated carbon determined from the analysis of pure-ethane adsorption at 273, 298, and 323 K. Figure 2. Fit to the pure-ethane experimental data of () Russell and LeVan 1997 using the PSD of Figure 1. Diamonds are for data at 273 K, open squares are for 298 K, and triangles are for 323 K. Ž. Russell and LeVan 1997 have reported data for the adsorption of pure water, pure ethane, and various mixtures of these components at temperatures of 273, 298 and 323 K. We have used the preceding method to calculate a PSD for BPL carbon based on their experimental data for ethane at all three temperatures. It was decided to use three temperatures, rather than a single temperature as in earlier work Ž. Davies and Seaton, 1999 , because it was intended to predict adsorption of waterrethane mixtures over this temperature range. After the single-pore isotherms for ethane were obtained from GCMC simulations, the AIE was solved using a Ž. nonlinear minimization routine MINOS, 1988 interfaced through the GAMS mathematical programming package Ž. GAMS, 1998 . The window of reliability was determined to lie between approximately 0.6 and 2 nm, and the maximum Ž resolution possible corresponding to the number of experi- . mental data points was used. The distribution with the optiŽ mal regularization the generalized cross-validation method and the L-curve yielded very similar regularization parame- . ters, resulting in essentially indistinguishable PSDs is shown in Figure 1, and the fit to the experimental data in Figure 2. The calculated PSD fits the experimental data well over the entire pressure range. Chemical Heterogeneity The PSD model for the structure of the carbon has been successful in predicting pureand multicomponent adsorpŽ. tion involving nonpolar Davies and Seaton, 1999 and weakly Ž. polar Heuchel et al., 1999 species. However, as we have discussed, a description of the chemical heterogeneity of the carbon is required when highly polar species, like water, are concerned. Here we attempt to include this extra dimension into our model carbon by placing oxygen-containing sites on the surface of the slit pores. The behavior of the real carbon is likely to be an average of the contributions of a wide variety of chemical groups Ž. Ž Boehm, 1994 . However, both experimental Barton et al., August 2003 Vol. 49, No. 8 AIChE Journal2062 .Ž. 1994 and molecular simulation Jorge et al., 2002 results suggest that the amount of water adsorbed on activated carbon depends primarily on the concentration of oxygen on the surface, regardless of its functionality. Therefore, we have assumed that all the polar sites present in the carbon can be represented by carbonyl groups. As for the location of these sites, it is known that they predominate near the edges of the graphite plates. Since our model represents the carbon structure by a collection of Ž. semi-infinite pores with no edges , we have chosen to place the polar sites on a square lattice superimposed on the pore walls. Furthermore, the sites were regularly distributed and separated by at least 0.75 nm. These simplifications will have an effect on the structure of the adsorbed water in the model, which will be different from that in real carbons. However, we have observed from previous molecular simulation studies that the distribution of polar sites on the surface of the carbon has a strong impact at low pressure, but almost no effect Ž. on the amount adsorbed at high pressures Jorge et al., 2002 . Ž An analysis of pure-water adsorption at low pressure prefer- . ably in the Henry’s law limit might, in principle, provide us with information on the site distribution. Unfortunately, such low-pressure experimental data are hard to come by and were unavailable for the particular carbon under study. However, since most of the available experimental data were obtained Ž. at high partial pressures of water see Figure 5 , the assumption of a regular site distribution in a pore with no edges should have a negligible effect on the results. Furthermore, the situation of technological interest, and difficulty, is when water adsorption on the carbon is substantial, which means that for most practical cases, the precise distribution of the polar sites on the carbon surface will not be an important factor. The last factor one must consider is the distribution of polar sites between pores of different width. An obvious first Ž choice is simply to say that this distribution is uniform that . is, all pores contain the same density of carbonyl groups . This assumption, together with the previous ones, reduces the process of characterizing the chemistry of the activated carbon to the determination of a single variableᎏthe concentraFigure 3. Pore-size distribution of BPL activated carbon reduced to 15 discrete pore sizes. tion of carbonyl groups on the surface. To estimate this concentration, we have first located the inflection point in the Ž. experimental typically S-shaped pure-water adsorption Ž. isotherm at 298 K see Figure 5 . Then we have determined the modal pore width of the PSD obtained from pure-ethane Ž. adsorption Figure 1 , which was calculated as 1.24 nm. A series of GCMC simulations in a model pore of that width, for different concentrations of polar sites, was performed, taking note of the variation of the pressure of spontaneous condensation with the site concentration. The chosen concentration was such that the pressure of spontaneous condensation in the single-pore isotherm corresponded roughly to the inflexion point of the experimental isotherm. This concentration, 2.214 ␮ molrm 2 , was used as a first estimate in the model. After choosing the concentration of carbonyl groups, we are able to calculate the total pure-water adsorption isotherm corresponding to the model carbon. This involves calculating a series of single-pore isotherms by GCMC, one for each pore size of the PSD, with the chosen concentration of polar sites. However, the effort involved in calculating 57 different water isotherms is prohibitive, even after fully optimizing the calculation of the long-ranged potential for charged molecules. We have therefore reduced the resolution of the original PSD to only 15 pores. This number was sufficient to retain all the characteristics of the PSD, while significantly reducing the computation cost. The reduced PSD is shown in Figure 3 for comparison. It is worth noting that the fit to the pure-ethane adsorption isotherms using this reduced PSD is essentially indistinguishable from that obtained with the original PSD Ž. Figure 2 . This is in agreement with the conclusions of Davies Ž. and Seaton 2000 , corroborating their statement that the resolution of the PSD is of minor importance when adsorption prediction is concerned. Some of the single-pore isotherms calculated from GCMC simulation are shown in Figure 4. Figure 4. GCMC adsorption isotherms for pure water at 298 K in slit-shaped pores with a concentration of surface sites of 2.214 ␮ molr rr rrm 2 and different widths. Filled squares ᎏ0.76 nm; open circlesᎏ0.84 nm; crosses ᎏ0.89 nm; open trianglesᎏ1.04 nm; filled circles ᎏ1.24 nm; open squaresᎏ1.44 nm; filled triangles ᎏ 2.04 nm. The pressure is reduced by the saturation pressure Ž. atthesametemperature 3167Pa. August 2003 Vol. 49, No. 8AIChE Journal 2063 We can see that at low pressure there is hardly any water adsorbed in the pores. At a certain characteristic pressure Ž. specific to each pore width , there is a steplike increase in the amount adsorbed, until the pores are completely filled. This discontinuous transition is evidence of capillary condensation, and was observed in all the isotherms at this temperature, down to the smallest pore width. It is worth noting that during a GCMC simulation in a semi-infinite pore at subcritical conditions, the true equilibrium vapor᎐liquid transition is never attained in practice. Instead, during adsorption calculations, the vaporlike phase remains metastable up to presŽ sures well above equilibrium liquid metastability also occurs .Ž . during desorption runs Jorge et al., 2002 . Thus, the vertical step in the single-pore isotherms represents a nonequilibrium w spontaneous condensation process. The equilibrium transiŽ tion can, however, be simulated using special methods Jorge .x and Seaton, 2002b . The vapor metastability that occurs in porous solids during real adsorption experiments is not strictly comparable to the corresponding phenomenon in a simulation, as the dynamics of the simulated and real processes are Ž. distinct Sarkisov and Monson, 2000; Schoen et al., 1989 . Ž. However, Sarkisov and Monson 2000 have observed that in a complex pore structure, the hysteresis arising from GCMC simulations seems to correspond well to the experimental hysteresis. Furthermore, recent studies in well-characterized porous materials involving mostly open unconnected pores Ž. Ravikovitch et al., 2001 show good agreement between simulation and experiment and suggest that the experimental adsorption branch is closer to the point of spontaneous condensation in GCMC simulation than to the equilibrium transition. In light of this, we will assume that the spontaneous condensation pressure obtained from GCMC is a good approximation to the value observed experimentally. The total adsorption isotherm for the model carbon was obtained by substituting the single-pore isotherms and the PSD into Eq. 1, and calculating the amount adsorbed, N. This isotherm is compared to the experimental data in Figure 5. As we can see from this figure, the simulated isotherm does not capture the S-shape of the experimental isotherm. Instead, it shows two distinct ‘‘steps,’’ reflecting the bimodal character of the PSD. Another important observation is that the simulation results lie almost entirely at lower pressures than the experiment, which suggests that the concentration of surface sites is too high. In a more encouraging note, the value of simulated adsorption at the bulk saturation pressure closely matches the experimental value, which lends further credibility to the representation of the structure of the carbon. We have repeated the procedure described earlier for a Ž 2 . lower concentration of carbonyl sites 1.476 ␮ molrm . The simulated isotherm is also shown in Figure 5. The reduction in the surface site concentration shifts the isotherm to higher pressures, but does not alter its shape. In fact, the simulation results will retain the bimodal character for any site concentration. Changing the type of oxygen-containing group used in the model, or changing the interaction parameters of the carbonyl group, would have a similar effect on the simulation results. Choosing a different distribution of the sites on the pore walls should bring about a significant change in the low-pressure region of the single-pore isotherms. However, the water uptake is mainly a result of capillary condensation, () Figure 5. Experimental data of Russell and LeVan 1997 () for pure water at 298 K diamonds and isotherms from the model carbon at two different concentrations of uniformly distributed 2 () carbonyl sites: 2.214 ␮ molr rr rrm solid line and 2 () 1.476 ␮ molr rr rrm dashed line . The dotted line connecting the experimental data points is a guide to the eye. which means that the total adsorption isotherm would be virtually unchanged. Since relaxing our two assumptions of uniform site type and regular site distribution on the pore walls will not improve the comparison between simulation and experiment, we are left with relaxing the third and lastᎏthe uniform site distribution on pores of different width. We have now allowed the density of carbonyl sites in each of the 15 pores to vary, while trying to match the simulation to the experimental Ž isotherm. This means using a site density distribution rather . than a single overall density value to represent the chemical heterogeneity of the carbon. Such a distribution was calculated by fitting the simulated isotherm to the experimental Ž pure-water adsorption data. The optimal distribution which . we will term distribution I is shown in Figure 6, and the corresponding fit to the experimental pure-water isotherm is shown in Figure 7. Combining a PSD for the structure of the carbon with a distribution of carbonyl sites to describe the chemical heteroŽ geneity provided good fits to both pure-ethane the fit shown in Figure 2 is virtually unaffected by the presence of polar . sites and pure-water adsorption isotherms. The distribution of polar sites shows two peaks, one around pores of 1.44 nm in width and the other around 0.76 nm. The latter peak is somewhat ill-defined, due to the shortage of experimental data at low pressure, and it is natural to ask how a distribution that omitted this peak altogether would perform. To find out, we have calculated isotherms for the model carbon using the two distributions of surface sites shown in Figure 8. Distribution II is based on a normal distribution centered around Ž 1.49 nm with a standard deviation of 0.3 nm thereby, incorporating the main peak of distribution I, but neglecting the . peak at low pore size , while distribution III is based on an Ž inverse exponential and is, therefore, monotonically increas- . ing, unlike both distributions I and II . When calculating these August 2003 Vol. 49, No. 8 AIChE Journal2064 Figure 6. Optimal distribution of polar sites on BPL carbon determined by analysis of pure-water ad- () sorption at 298 K distribution I . distributions, the density of surface sites was rounded off to allow for an integer number of carbonyl groups in the simulation cell. In Figure 9 we compare the experimental isotherm with the simulated one, obtained using distributions II and III. As we can see, distribution II fits the experimental data very satisfactorily above a relative pressure of 0.4. This is expected, since the peak of distribution II closely resembles the second Ž. peak of distribution I the best-fit distribution . This peak lies in the region of large micropores, which are responsible for the water uptake at high relative pressures. However, the low-pressure adsorption is severely underestimated; this can also be related to the polar site distribution. The low-pressure uptake is mainly dictated by adsorption in the small micropores. In distribution II, these pores are almost devoid of polar sites, and, thus, condensation occurs much later. The fit to experiment is much worse for the case of distribution Figure 7. Fit to the pure-water experimental data of () Russell and LeVan 1997 at 298 K on BPL carbon using the PSD of Figure 1 and the polar site distribution of Figure 6. The diamonds represent the experimental data and the solid line is the simulated isotherm. () Figure 8. Two ‘‘smooth’’polar-site distributions: a dis- () tribution II; b distribution III. III. In fact, it is almost as bad as for the uniform site distributions shown previously. The reason for this is that, because it is monotonically increasing, distribution III is unable to com- () Figure 9. Experimental data of Russell and LeVan 1997 () for pure water at 298 K diamonds and isotherms from the model carbon with two different distributions of carbonyl sites: distribu- ()( ) tion II solid line and III dashed line . August 2003 Vol. 49, No. 8AIChE Journal 2065 pensate for the bimodal PSD, with the simulated isotherm retaining the two distinct steps. From this analysis, we can conclude that while a suitably chosen unimodal distribution of surface sites can describe the data above a relative pressure of around 0.4, below this pressure the peak at small pore size shown in distribution I must be taken into account. A precise description of the site density in these small pores is difficult in the present case, due to the reduced amount of experimental data obtained at low relative pressure. In summary, the predicted adsorption of water is quite sensitive to the precise form of the polar site distribution as a function of pore size. We are now able to compare the physical properties of our Ž model carbon composed of the PSD of Figure 3 and the . polar site distribution of Figure 6 with those of real BPL carbon, determined experimentally. We have calculated the Ž. micropore volume, the specific surface area Sand the surŽ. face oxygen concentration Cfor the model carbon. The O total volume of micropores for the model carbon is given by m Vsfw ␦ w3 Ž. Ž. Ý mic nn ns1 In order to compare this volume to the experimentally determined values on the same basis, the result of Eq. 3 was modified to account for the volume of the pore that is inaccessible wŽ. x to the adsorbate molecules see Jorge 2003 , for details . The specific surface area for the model is simply given by m2fw ␦ w Ž. nn Ss4 Ž. Ýw n ns1 Finally, using the polar site distribution and taking into account that each carbonyl group contains one oxygen atom, the total concentration of oxygen on the model carbon can be obtained from m2fw ␦ w Ž. nn Cscw 5 Ž. Ž. Ý O n w n ns1 Ž. where cw is the molar concentration of carbonyl sites per unit surface area in a pore of width w. Table 2 shows a comparison between the specific properties of the model with those of real BPL carbon. Bradley and Ž. Rand 1995 have obtained a consistent value for the total micropore volume of BPL carbon by applying the DubininRadushkevitch equation to isotherms of several nonpolar adTable 2. Specific Properties of BPL Carbon from Experimental Measurements vs. Those for the Model Carbon Model Experiment U 3 Ž. Vcm rg 0.417 0.43 mic UU 2 Ž. Smrg 925 1,088y1,120 † Ž. Cmmolrg 1.16 1.14y1.40 O U Ž. Taken from Bradley and Rand 1995 . UU Ž. Ž. Taken from Barton et al. 1998 , MacDonald et al. 2000 , and MacŽ. Donald and Evans 2002 . † Ž. Ž. Taken from Barton et al. 1996, 1997 , MacDonald et al. 2000 , and Ž. MacDonald and Evans 2002 . sorbates. Barton and coworkers have performed extensive experimental characterization studies on BPL carbon. They have obtained specific surface areas by applying the standard Ž BET method to N adsorption Barton et al., 1998; MacDon2 . ald et al., 2000; MacDonald and Evans, 2002 . They have also performed a series of temperature-programmed desorption studies to determine the oxygen concentration on the surface Ž of BPL carbon Barton et al., 1996; Barton et al., 1997; Mac- . Donald et al., 2000; MacDonald and Evans, 2002 . In Table 2, the range of surface areas and oxygen concentrations obtained on different samples is given. The micropore volume of the model carbon shows good agreement with the experimental determination, but the specific surface area is slightly below the BET surface area. However, it is likely that the BET method overestimates the real specific area of the carbon, since it does not take into account the enhancement of the adsorption energy in small pores. As for the oxygen concentration of the model, it is within the range of independently determined values. This comparison should only be taken as a rough guide, since the oxygen atoms included in the model effectively account for all the different sites of a real carbon that interact with water, including other heteroatoms and basic sites. Nevertheless, the good agreement obtained suggests that the approach presented here is physically reasonable. Therefore, it can be concluded that both the structure and the chemistry of the model carbon are consistent with experimental data obtained on real BPL carbon. Prediction of Multicomponent Adsorption We have assessed the predictive capabilities of the model by comparing simulation results with experimental measurements of a mixture of water and ethane on the same sample of BPL. The adsorption isotherm for the model was obtained by performing GCMC simulations in the 15 slit-shaped pores considered previously, with a concentration of carbonyl sites given by the distribution of Figure 6. These isotherms were then weighted by the PSD of Figure 3, according to Eq. 2, to yield the total adsorption isotherm. Figure 10 shows a comparison between the multicomponent adsorption isotherms obtained from experiment and from molecular simulation. The agreement is reasonable, even though there is a slight underestimation of water uptake at Ž high partial pressures of ethane corresponding to low mole . fractions of water in the bulk mixture . This is balanced by an overestimation of the ethane uptake in this region. Water uptake at very low mole fractions of water is mainly dictated by adsorption in pores of width below 0.84 nm, since ethane is preferentially adsorbed in the larger pores. This is precisely the range of pore sizes in which the polar site distribution is poorly defined. As explained previously, this is because of the restricted amount of available experimental pure-water data at low relative pressure. It is, therefore, not surprising that the multicomponent adsorption predictions in this range, where the polar site density has a particularly strong effect on adsorption, are poorer. The different adsorption behavior in pores of different sizes is illustrated in Figure 11, which shows snapshots obtained from waterrethane simulations in two pores of different width and at a high ethane partial August 2003 Vol. 49, No. 8 AIChE Journal2066 Figure 10. Experimental data of Russell and LeVan () ( ) 1997 for binary water solid diamonds and () ethane open diamonds adsorption at 298 K ( with isotherms from the model carbon solid ) line using the PSD of Figure 3 and the polar-site distribution of Figure 6. pressure. As we can see, the smaller pore is filled with water, while the large pore contains mostly ethane molecules. Again, it would be interesting to see what the behavior of the model carbon would be, with respect to waterrethane mixture adsorption, if a uniform carbonyl site distribution was Ž. assumed instead of distribution I . We have calculated the Figure 11. Molecular simulation of a mixture of ethane () partial pressure s73,416 Pa and water () partial pressure s1584 Pa at 298 K adsorbed on activated carbon pores of width: () () a 0.8 nm; b 1.24 nm. The magenta spheres represent the L-J centers of ethane, the transparent gray spheres represent the L-J centers of water, and the green spheres represent the carbonyl sites. The point charges of the models are represented by solid spheres, blue for a negative charge and red for positive. Figure 12. Experimental data of Russell and LeVan () ( ) 1997 for binary water solid diamonds and () ethane open diamonds adsorption at 298 K with isotherms from the model carbon using the PSD of Figure 3 and two different concentrations of uniformly distributed carbonyl 2 () sites: 2.214 ␮ molr rr rrm solid lines and 1.476 2 () ␮ molr rr rrm dashed lines . multicomponent adsorption isotherms at the same conditions, using two different values for the uniform surface site density in the pores. The results are compared with the experiment in Figure 12. As we can see, the predictions are now extremely poor, with water adsorption being substantially overestimated and ethane adsorption underestimated. This is due to the large site density in the pores around 1 nm, which leads to a high water uptake in the entire pressure range. The effect here is analogous to the overestimation of the pure-water adsorption at relative pressures below 0.5, visible in Figure 5. From this analysis, the importance of a complete and detailed description of the chemical heterogeneity of the activated carbon is once again emphasized. We have also assessed the performance of two classical thermodynamic methods for this system, namely the IAST Ž. Myers and Prausnitz, 1965 and the Doong-Yang method Ž. 1988 . The implementation of these methods was as described in the original articles, except where indicated below. In the case of the IAST, the integration of the pure-component isotherms was performed numerically by Simpson’s rule. When extrapolation of the pure-ethane isotherm to high pressure was required, this was done using the Toth equation Ž. Toth, 1962 . The nonideality of the bulk gas phase was taken into account by calculating fugacities using the Peng-Robinson equation of state. In the case of the Doong-Yang method, the Dubinin-Astakhov equation was fitted to the pure-component isotherms. Although the method allows for different limiting pore volumes for each of the components, we have Ž used the same value for both water and ethane using differ- . ent volumes yielded no improvement to the predictions . The predictions using IAST and the Doong-Yang method are presented in Figure 13. As we can see, both methods fail to correctly predict the binary water adsorption isotherm, substantially overestimating the uptake. The IAST substantially overpredicts the adAugust 2003 Vol. 49, No. 8AIChE Journal 2067