Adsorption on heterogeneous surfaces. Calculation of the adsorption

Adsorption of Spherical Molecules in Probing the Surface Topography: 2. ... Numerical Method for Analysis of Surface Heterogeneity in a Case of Finite...
1 downloads 0 Views 1MB Size
Download PDF
Langmuir 1993,9, 2593-2605

2593

Adsorption on Heterogeneous Surfaces. Calculation of the Adsorption Energy Distribution Function or the Affinity Spectrum Luuk K. Koopal' and Cornelis H.W.Vos Department of Physical and Colloid Chemistry, Wageningen Agricultural University, Dreyenplein 6,6703HB Wageningen, The Netherlands Received August 3,1992. I n Final Form: November 30, 199.P A brief overview is given on the calculation of adsorption energy distribution functions or affinity spectra on the basis of adsorption or binding data. First the local isotherm approximation or LIA methods are discussed. These methods lead to expressions for the distribution function in terms of derivatives of the overall isotherm which makes smoothing of the isotherm data of particular importance. Good results are obtained with a smoothing spline technique in combination with generalized cross validation and some physical constraints to determine the smoothingparameter. In the second part of the paper three numerical methods, viz., CAESAR, REMEDI, and CAEDMON-W, are compared. CAESAR is an existing program, R E m D I is essentially similar to the regularization method proposed by House, and CAEDMON-W is a modification of the program CAEDMON of Sacher and Morrison. The selection of the parameter which stabilizes the solution is made in a similar way in all methods, to be able to make a precise comparison between the methods. Also the same local isotherm equations are used. The methods are tested with simulated isotherm data, and the effects of random and systematic errors are investigated. Application of the methods to Nz adsorption data on a series of carbon blacks and a series of modified silicas shows that the results obtained with the different methods are very similar. It may be concluded that the methods work very well and that spurious oscillations are suppressed. Some preference is to be given to CAESAR or REMEDI, because with CAEDMON-W only a limited number of points of the distribution function can be calculated.

Introduction Adsorption by heterogeneous surfaces and/or binding of small molecules to heterogeneous ligands with a multiplicity of specific binding sites is important in many fields of science, for instance, in the field of chemical technology, the adsorption of trace components to activated carbons and the characterization of adsorbent surfaces,l in biochemistry, the binding of steroids to tumors,2 in immunochemistry, the binding of antibodies to heterogeneous multivalent antigen systems,314 and in environmental sciences, the binding of trace amounts of toxic components to natural particles such as metal hydroxides and humic It is generally recognized that the heterogeneity of substrates or adsorbents with respect to binding is rather important, but very little is known about the heterogeneity. Therefore, it is important to have methods available for the analysis of binding or adsorption data that lead to information on the heterogeneity or, more precisely, on the distribution of the affinity constants or the standard state Gibbs energies of adsorption. In biochemical, environmental, and immunological literature the distribution of affinity constants is often called "the affinity spectrum";' in the adsorption literature, one generally refers to it as "the adsorption energy di~tribution".~*Q Gas adsorption has for a long time been a standard Abstract published in Advance ACSAbstracts, Auguat 15,1993. (1) McEnaney, B.; Mays, T. J.; Cawton, P. D. Langmuir 1987,3,695. (2) Mechanick, J. I.; Peakin, C. S. Anal. Biochem. 1986,157, 221. (3) Gandolfi, A.; Strom, R.J. Theor. Biol. 1981, 92,57. (4) Sciutto,E.; Garat, B.;Ortega, E.; Larralde, C.Mol. Zmmunol. 1987, 24,577. (6) Gamble, D. S.; Langford, C. H. Enuiron. Sci. Technol. 1988, 22, 1925.

(6) Buffle, J.; Altmann, R. S.;Filella, M.; Tessier, A. Geochim. Cosmochim. Acta 1990,64, 1535. (7) Thakur,A. K.; Munson, P. J.; Hunston, D. J.; Robard, D. Anal. Biochem. 1980,103, 240. (8) House, W. A. Colloid Science. SpecialistPeriodicalReport;Everett, D. H., Ed.; Royal Society of Chemistry: London, 1983; Vol. 4, Chapter 1.

technique to study solid adsorbents and to characterize surfaces with respect to their area and porosity (see e.g. ref 10). A relatively high degree of sophistication has been reached with the measurement of nitrogen, krypton, and argon adsorption. Commercialequipment is available that allows accurate measurements in the low pressure range, or the (sub)monolayer domain of adsorption. Data in this range are most suited for analysis of heterogeneity. In the field of biochemistry and environmental chemistry, when information is sought about specific substrate ligand complexes,it is often difficult to obtain the complete binding curve, i.e., from an entirely unoccupied substrate to a substrate on which all binding sites are occupied, and only a "window" of the full range can be measured. Moreover, it may not be possible to obtain highly accurate data. Nevertheless, the need for information on the affinity distribution is still present, and methods are required that can deal with this problem. The fact that different types and qualities of data sets have to be analyzed has resulted in different methods for the calculation of the affinity distribution function being proposed. In the present paper we will first describe some general aspects of the heterogeneity analysis. A brief review of the existing methods to analyze adsorption isotherms or binding curves for heterogeneity will be given with special attention to the semianalytical methods. These methods are especially suited when only a window of the full data range is present. In the second part of the paper three sophisticated numerical methods developed for the analysis of gas adsorption isotherms will be compared. The terminology from the adsorption literature will be used, but the ideas relate equally well to the binding of small molecules to polyfunctional ligands. (9) Jaroniec, M.; Madey, R. Physical Adsorption on Heterogeneous Solids; Elsevier: Amsterdam, 1988. (10) Gregg,S . J.; Sing,K. S. W .Adsorption Surface Area and Porosity, 2nd ed.; Academic Preae: London, 1982.

0743-7463/93/2409-2593$04.00/00 1993 American Chemical Society

Koopal and Vos

2594 Langmuir, Vol. 9, No. 10, 1993

General Aspects of the Heterogeneity Analysis Adsorption on a heterogeneous surface can be considered as binding to a collection of various homogeneous subsurfaces, the homogeneous subsurface being an ensemble of sites for which the standard Gibbs energy of adsorption, or the intrinsic affinity, is a constant. For a continuous distribution of affinities the summation of the adsorption over all homogeneous subsurfaces can be done by integrating over all possible values of the affinity, leading to e,(x) = J’e(K,c,x) m g K ) doog K )

(la)

In eq l a e&) is the normalized overall adsorption as a function of the activity x at a given temperature, B(K,C,x) is the “local” isotherm describing the adsorption on an ensemble of “equal intrinsic affinity” sites, f(1og K ) is the distribution of log K, and A is the range of the log Kvalues. The local isotherm is a function of the intrinsic affinity K, the lateral interactions or “nonideality” expressed by C, and the bulk activity x of the adsorbing component. One may also write eq l a in terms of free energies by realizing that log K = Agao/RT:

..

= J&u,c,x)f(wd U

Ob)

where Agao/RTis abbreviated as U. Note that in principle Agao/RT includes both an energy and an entropy contribution. In the case of gas adsorption, the thermal entropy contribution is small and often neglected. In the case of adsorption from solution, the nonconfigurational entropy can be substantial (e.g., due to hydrophobic bonding) and &a0/RT is a true free energy. Only for the sake of convenience we have abbreviated Agao/RT as U,but this does not imply that for adsorption from solution the nonconfigurational entropy contribution can be neglected. In eq 1 the adsorption is expressed as a dimensionless quantity which can be obtained from the measured adsorption by dividing this quantity by a normalization constant (usually the maximum adsorption in a monolayer). It should be realized that eq 1 is based on a physical model of the adsorption in which it is assumed that C and K vary independently. This is equivalent to assuming that the local isotherm and the distribution function are entirely independent. The exact nature of the forces which lead to the distribution of log K or U is irrelevant. The distribution might be due to the presence of surface sites which differ in chemistry, but it might equally well originate from a difference in local structure, such as the presence of micropores of different sizes, or be due to a combination of factors. Knowledge of these factors and the interaction mechanisms may be helpful in establishing an appropriate expression for 8or in interpreting the result once a distribution function has been calculated. When e&) is measured, eq 1can, in principle, be solved for the distribution function provided an expression is available for the local isotherm 8. The expression for the local isotherm determines which interactions are included in log K or in U. In the presence of lateral interactions the standard Gibbs energy of adsorption or intrinsic affinity, representing the surface-adsorbate interactions, is not the only contribution to the adsorption Gibbs energy or total affinity. The lateral interactions between the adsorbate molecules (indicated by C) contribute also and have to be accounted for separately, if possible. The way this is done depends on the type of heterogeneity. General practice is to consider two limiting situations: random and patchwise heter~geneity.~?~ For random heterogeneity the lateral interactions are averaged over the entire surface.

For patchwise heterogeneity the lateral interactions are averaged per patch only. The patches are assumed to be large enough to neglect boundary effects. By incorporating the adsorbate-adsorbate interactions separately from the adsorbate-adsorbent interactions, it is possible to retain the standard Gibbs energy of adsorption (the intrinsic affinity) for each homogeneous subsurface as a separate parameter. When lateral interactions are not considered explicitly, an “apparent” affinity distribution will be found. In the biochemical and environmental literature, where rather complicated substrates and/or ligands are present, it is very difficult to fully separate ligand-ligand and ligandsubstrate interactions. In such cases one has to be satisfied with an apparent affmity distribution only. With respect to the technique of the heterogeneity analysis the separation between adsorbate-adsorbate and adsorbate-surface interactions is immaterial (provided C andKare independent). Once agreementhas been reached about the equation for the local isotherm, eq 1 can be inverted and an affinity distribution can be found. The advantage of separating adsorbate-adsorbate interactions from adsorbate-surface interactions lies in the fact that an intrinsic affinity distribution provides direct information on the solid-adsorbate interactions. For instance, studying metal ion adsorption to a humic acid on the basis of an intrinsic affinity distribution allows coupling of the observed log K values to log K values (stability constants) for low molecular weight complexes and, hence, to the chemistry of the reactive groups. Apparent affinities do not discriminate between adsorbate-surface interactions and adsorbate-adsorbate interactions. Therefore, it is more difficult to understand their meaning. An apparent affinity distribution is essentially an alternative way of representing the adsorption or binding isotherm. Whenever the affinity representation allows more insight in the binding behavior than the isotherm itself, it will be useful to calculate the apparent distribution function. Both the calculated intrinsic and apparent affinity distributions are useful aids for the selection of an adequate equation for the overall adsorption isotherm. A general expression for the local isotherm is the following equation: (2) e(K,c,x) = ~ c x / (+i K C ~ ) For C = 1 this equation is the well-known Langmuir equation. In general however C # 1. Expressions for the function C can be derived from adsorption models for homogeneous surfaces.ll For instance, a mean field approximation taking into account nearest-neighbor interactions only leads to the Fnunkin-Fowler-Guggenheim (FFG) equation12J3and

C = exp(-wB*/RT) (3) where w represents the net standard Gibbs energy of lateral interaction and R and T have their usual meanings. For random heterogeneity 8* = et; for patchwiseheterogeneity 8* = 8. Preferably w should be obtained independently from other experiments and/or model considerations; w depends mainly on the adsorbate, and it may contain both enthalpic and entropic contributions. For charged surfaces C also includes a Boltzmann factor based on the surface potential $8 to account for the Coulombic interactions.14 (11) Rose, S.;Olivier, J. P. On Physical Adsorption; Wiley-Interscience: New York, 1964. (12) Fmkin, A. N. 2.Phys. Chem. 1926,116,466. (13) Fowler, R.H.;Guggenheim, E, A. Statistical Thermodynamics; Cambridge University Press: Cambridge, 1966; p 426.

Adsorption on Heterogeneous Surfaces

In the case of mobile adsorption, a Hill-de Boer-type equation should resultl5and C includes not only the lateral interaction expression but also a factor accounting for the difference in configuration entropy between localized and mobile adsorption. A virial-type local isotherm16results if eq 3 is expanded in a series, the virial coefficients depending on the expression for C . For random heterogeneity, C is a function of Bt only and, provided C is known, the overall adsorption can also be expressed as a function of Cx instead of x . By doing this, the local isotherm should be expressed as 8(K,Cx) which transforms eq 2 to the simple Langmuir equation. For complex heterogeneous ligands occurring in natural systems it is often a good first-order approximation to assume that the surface is random heterogeneous, so that this type of transformation of the data can be made. A check on the assumption is to plot overall isotherms measured under different conditions (different values of C ) as a function of Cx. Merging of thus obtained isotherms in a master curve indicates the suitability of the model used to calculate C ; see for instance refs 17 and 18. For patchwise heterogeneity such a simple transformation is not possible because C is in this case a function of 8. In order to solve eq 1 for f(1og K ) or f(i3for known values of et, it should be realized that eq 1is a Fredholm equation of the first kind. Solution of the Fredholm equation with respect to f(1og K ) when both and 8(K,C,x) are known is not a trivial task because a small variation in the data 8&) can cause a large variation in f(1og K).8Jg Since data are subject to experimental error, we are faced with the possibility of finding a solution whose physical significance is obscured by spurious oscillations due to experimental error. In the literature two main routes have been followed to solve eq 1for f(logK),without making assumptions about the distribution function. In the first group of methods an approximationis made which makes it possible to solve eq 1 analytically for the distribution function. A more rigorous alternative is to use a special numerical method, in which the ill-posedness of the Fredholm equation is addressed explicitly. In the next section the semianalyticalmethod will be briefly reviewed;thereafter we will return to the purely numerical methods.

Local Isotherm Approximation (LIA) Methods In view of obtaining an anlytical expression for the distribution function several researchers have made approximations of essentially the local isotherm function. Examples of these methods are the condensation approximation (CA), which represents the isotherm by a step function,20121the asymptotically correct approximation (ACA),based on a linear i~otherm,2~9~3and the logarithmic symmetrical approximation (LOGA), based on a function which can approximate most local isotherm equation^.^?^^ LOGA's mostly used form is the LOGAl which describes (14)Koopal, L. K. In Colloid Chemistry in Mineral Processing; Laskowski, J. S., Ralston, J., Eds.; Elsevier: Amsterdam, 1992; Chapter 2, p 37. (15) Hill, T. L. Introduction to Statistical Thermodynamics, 2nd printing; Addison-Wesley: Reading, MA, 1962. (16) Morrbon, I. D.; Ross, 5.Surf. Sci. 1973, 39, 21. (17)De Wit, J. C. M.; Van Riemadijk, W. H.; Nederlof, M. M.; Kmniburgh, D. G.;Koopal, L. K. Anal. Chim. Acta 1990,232, 189. (18) De Wit, J. C.M.;VanRiemsdijk, W.H.;Koopal,L.K. Water,Air, Sod Poliut. 1991,57-58, 339. (19) Men, P. H. J. Comput. Phys. 1980,38,64. (20) Roginski, S. S. Adsorption and Catalysis on Heterogeneow Surfaces; Academy of Sciencea of USSR Moscow, 1948. (21) Harris, L. B. Surf. Sci. 1968, 10, 129. (22) Cerofolini, G.F. Surf. Sci. 1971,24, 391. (23) Cerofolini, G. F. Thin Solid Films 1974,23, 129.

Langmuir, Vol. 9, No. 10, 1993 2595

a Langmuir-type function quite closely. All LIA methods lead to expressions for the distribution function in terms of derivatives of the overall isotherm. Nederlof et a1.W have shown that analyticalexpressions for the distribution function taking the form of a series of derivatives of 8,(x) can be considered as solutions based on a local isotherm approximation. Examples of methods which implicitly belong to the LIA methods are the Rudzinski-Jagiello (RJ)method,26*nthe first- and secondorder affinity spectrum methods (AS1 and AS2),7*aand the differentialequilibrium function (DEF)method.SaThe RJ method is well known in the adsorption literature; the AS methods and the DEF method stem from the biochemical and environmental literature. The AS methods essentially go back to the work of Schwarzl and StavermanmIsoand Ninomiya and Ferry.31 The DEF method is introduced by Gamble and co-worker~S*~~ and extended by Buffle and his gr0up.~9~~9" Of these methods the RJ and the AS2 methods turn out to be specific cases of the LOGAl method, whereas the AS1 method closely resembles the CA method.24125 The DEF method has some features in common with the ACA method, but the second part of the local isotherm corresponding to the DEF method is physically impossible.25 Although the LIA and related methods lead to an analytical expression for the distribution function, part of the instability problem still exists. This is due to the fact that the distribution is obtained by taking the first and higher derivatives of the overall adsorption isotherm. Taking derivatives of experimental data without preprocessing the data easily leads to spurious peaks resulting from experimental errors. Thus, in order to apply the semianalyticalmethods, a data processing routine should be available which deals with the error in experimental data. Nederlof et ale3" have discussed this problem in more detail and have tackled it by using a smoothing spline technique. For known error the variance in the data can be used to find a proper smoothing parameter. In addition to the variance, an expression is required for the goodness of fit and for the smoothness of the spline as a function of the smoothing parameter. For unknown error, a situation often occurring in practice, a suitable criterion to determine the smoothing parameter is provided by the generalized crow validation (GCV) technique38939 in combination with physical constraints that must be obeyed. Once a proper smoothing parameter is determined, it is possible to represent the data with a spline and to calculate (24) Nederlof, M. M.; Van Riemsdijk, W. H.; Koopal, L. K. Environ. Sci. Technol. 1992, 26, 763. (25) Nederlof. M. M.: Van Riemsdiik, W. H.:. K w -~ a.lL. . K. J. Colloid Interface Sci. 1990, 136,410. (26) Rudzmki, W.; Jagiello, J.; Grillet, Y.J. Colloid Interface Sci. 1982, 87, 478. (27) Jagiello, J.; Ligner, G.;Papirer, E. J. Colloid Interface Sci. 1989, 137, 128. (28) Hunston, D. L. Anal. Biochem. 1975,63,99. (29) Schwarzl, F.; Staverman, A. J. Physica 1952,18,791. (30) Schwanl, F.; Staverman, A. J. Appl. Sci. Res. 1963,4, 127. (31) Nmomiya, K.; Ferry, J. D. J. Colloid Interface Sci. 1969,14,36. (32) Gamble, D. S. Can. J. Chem. 1970,48, 2662. (33) Buffle, J. Complexation Reactioha in Aquatic Systems: An Analytical Approach; Ellis Horwood: Chichester, U.K., 1988. (34) Altmann, R. S.; Buffle, J. Geochim. Cosmochim. Acta 1988,52, 1505. (35) Nederlof, M. M.; Van Riemsdijk, W. H.; Koopal, L. K. In Heavy

Metals in the Environment; Vernet, J.-P., Ed.; Elsevier: Amsterdam, 1991; p 365. (36) Nederlof, M. M.; Van Riemsdijk, W. H.; Koopal, L. K. Submitted for publication. (37) Nederlof, M. M. Analyaisof BindingHeterogeneity. Ph.D.Theaie, WageningenAgriculturalUnivereity, Wageningen,The Netherlands, 1992. (38) Silverman, B. W. J. R. Stat. SOC.B 1986,47,1. (39) Craven, P.; Wahba, G.Numer. Math. 1979, 31, 377.

2596 Langmuir, VoE. 9, No. 10, 1993

Koopal and Vos

the distribution function on the basis of that spline. In addition, from the variance of the spline, the uncertainty in the distribution function can be indicated by confidence intervals. The latter indicate whether peaks found are significant. In general GCV alone is not enough to obtain a proper smoothing parameter. This has to do with the fact that the error in the data is often not entirely random and/or that the data set is of limited size. Additional information is provided by the physical constraints. These constraints can be derived from the properties of both the overall and the localisotherms. Since, according to the local isotherm, adsorption on a homogeneous patch increases with increasing x , it follows that the overall isotherm should also increase with x . Therefore, a generally valid constraint is that the first derivative of the measured isotherm should be positive for all values of x . Higher order constraints depend on the local isotherm function used. Apart from the fact that confidence intervals can be calculated, a main advantage of the procedure described by Nederlof et al.3M7 is that the smoothingparameter is obtained in an objective way. For an optimum application of the GCV criterion it is essential to perform the smoothing procedure directly on the signal that is measured and not on an already "transformed" data set, unless the propagation of error is properly accounted for by weighting factors in the spline analysis. Examples of analyses of isotherms can be found in the work of Nederlof et al.3M7 Data with a rather large error will always lead to rather smooth distribution functions due to the necessary smoothing. Small peaks in a distribution are hard to recover with the LIA methods when the peaks are present on a "background" of a continuous distribution. A small peak at the high affinity tail of a distribution can be discovered fairly well. For systemsfor which the local isotherm equation is not known precisely (which is the case for many systems) the LIA methods are well suited to calculate the distribution function. This is especially the case when only a window of data is present. When general agreement exists about the precise form of the local isotherm, the LIA methods will always show a lower resolution than an adequate numerical method.40 A main advantage of the LIA methods is that they show the relation between the overall adsorption isotherm and the affinity distribution function.

Numerical Methods The second class of methods to solve the integral equation without making a priori assumptionswith respect to the distribution are numerical techniques. In principle these methods can be used with a variety of local isotherm equations. The first numerical technique used in the field of adsorption science has been developed by House and Jaycock41~42and is called HILDAT2(heterogeneity investigated by a Loughborough distribution analysis). HILDA is a numerical refinement of a graphical iteration method proposed by Adamson and Ling.43 A second method with the acronym CAEDMON (computed energy distribution the monolayer) was initially developed by Ross and Morrison44and improved by Sacher ~~

~~

~

~

~

(40)Koopal,L. K.;Nederlof, M. M.; Van Riemsdijk,W. H. h o g . Colloid Polym. Sci. 1990, 82, 19. (41) House, W. A.; Jaycock, M. J. J. Chem. SOC.,Faraday Trans. 1 1977, 73, 942. (42) House, W. A.; Jaycock, M. J. Colloid Polym. Sci. 1978,256, 52. (43) Adameon, A. W.; Ling, I. Adu. Chem. Ser. 1961,33,51. See also: Adamson, A. W. Physical Chemistry of Surfaces, 4th ed.; Interscience: New York, 1982.

and Morrison.& Both in HILDA and CAEDMON the ill-posed-ness of the problem has not been addressed explicitly. In the adsorption literature House was the first to address the instability and applied a regularization technique proposed by PhillipsY and Twomey49 (PTH regularization). In principle in the regularization methods an extra condition is made to the distribution function f(0, namely, that f(0must be smooth. The extent of smoothing is governed by the order of the regularization and the regularization parameter. The PTH regularization is second order. Merz,19 Brown and cow o r k e r ~ ,and ~ * Butler ~ ~ et aLS2studied ways to determine the regularization parameter "automatically". Another technique, similar to regularization, which also accounts for the ill-posed-ness of the problem is singular value decomposition (SVD). In the immunological literature Hansonm applied this technique to hapten binding. Koopal and V0sM96Sdeveloped a method named CAESAR (computed adsorptive energies, SVD analysis result) to analyze gas adsorption results. The extent of smoothing in CAESAR depends on the errors present in the data and those introduced by the chosen local isotherm equation which, in general, is only approximately correct. Recently Von Szombathely et al.58 have critically reviewed the regularization methods (including SVD)and examined the general principles of regularization for solving the present problem. A program, INTEG, based on regularization, was written for the calculation of the distribution function. Widely used for similar problems in the field of physical sciences is the program CONTIN of Proven~her.'~J~J~ The CONTIN method is also based on regularization and can be used for inverting noisy linear operator equations in general. CONTIN has not yet been applied to adsorption problems. An entirely different method to solve f(0is to use the Fourier transform (FT) method. This method has been applied by Metcalfe,67 by Mechanick and Peskin? and most elaborately by Lum Wang and White.a Also with the FT technique problems will occur because of the error in the data. Lum Wang and White used a frequency extrapolation technique to minimize the effect of noise on the convoluted distribution function. The FT method is especially suited for very large data sets. In the literature, results obtained with the abovementioned methods have been compared, but due to differences in the choice of the local isotherm and the boundary conditions these comparisons could only be (44)Ross, S.; Morrison, I. D. Surf. Sci. 1975,52, 103. (45) Sacher, R. S.; Morrison, I. D. J.Colloid Interface Sci. 1979, 70, 153. (46) House, W. A. J. Colloid Interface Sci. 1978, 67, 166. (47) House, W. A. J. Chem. Soc., Faraday Trons 1 1978, 74, 1045. (48)Phillips, D. L. J. Assoc. Comput. Mach. 1962,9,84. (49) "homey, S. J. A8soc. Comput. Mach. 1963,10,97. (50) Brown, L. F.;Travis,B. J. InFundamentalsof Adsorption; Myers, A. L., Belfort, G., Eds.; Engineering Foundation: New York, 1984; p 125. (51) Britten, J. A.; Travis, B. J.; Brown, L. F. AZChESymp. Ser. 1984, 79, 7.

(52) Butler, J. P.; Rees, J. A.; Dawson, S. V. SZAM J. Numer. Anal. 1981, 18, 381. (53) Hanson, R. J. SIAM J. Numer. Anal. 1971,8, 616. (54) Koopal, L. K.; Vos, C. H. W. Colloids Surf. 1986,14, 87. (55) VOS,C. H. W.; Koopal, L. K. J. Colloid Interface Sci. 1986,105, 183. (56) Von Szombathely, M.;Bratier, P.;Jaroniec, M. J.Comput. Chem. 1992, 13, 17. (57) Metcalfe, I. M. Computational Advances in Colloid and Surface Science, Part 111. Ph.D. Thesis, University of Melbourne, Melbourne, Australia, 1985. (58) Lum Wan,J. A.; White,L. R. J.Chem. Soc., Faraday Trans. 1991, 87, 3051.

Langmuir, Vol. 9, No. 10, 1993 2597

Adsorption on Heterogeneous Surfaces

superficial. In the next section a rigorous comparison will be made between CAESAR and somewhat modified versions of the PTH regularization and CAEDMON. The modifications were made to enable a more rigorous comparison of the methods. Essentially the basic algorithm of each method is used, but the local isotherm equation and the boundary conditions are adapted to those used in CAESAR. Moreover, f (v) is only allowed to take non-negative values. With these modifications it is possible to use the same criterion for all methods to select a value for the parameter which controls the spurious oscillations in the distribution function. This aspect and the use of the same local isotherm make a rigorous numerical comparison of the three methods possible. For numerical evaluation eq 1is in general approximated by a summation: n

non-negative least-squares (NNLS) procedure of Lawson and Hansons9 which is especially suited for the present problem.4s The numerical rank r of A is, in practice, the number of rows or columns of matrix A which may be considered as linearly independent in the numerical analysis. The value of r will depend on the machine precision, the experimental error in et, and the approximation that is inherent to the choice of a specific local isotherm function. In practice the machine precision determines the lower limit of qr (10-6 in CAESAR). Starting with this value, the vector f, is determined. Then f, is multiplied by the original matrix A to obtain et,, = Afc (9) To check this solution, the deviations between & ( x j ) and &,,(Xj) have to be considered. The solution f, is required such that m

where n (Im)is the chosen number of subintervals on the energy interval, A, each of equal width AU = ( U m h- Urns)/ n, m is the total number of experimental points of the isotherm (restricted to the monolayer region), e(Ui,xj) is the value of e( U,X)a t U = Ui and x = X j , f (Vi) is the value is the value of e&) a t x = X j . off( v) at U Vi,and N is a normalization constant equal to 1when the overall isotherm is expressed as a fractional coverage, or equal to the maximum monolayer adsorption when the overall adsorption is given in adsorption units. The set of linear equations derived from eq 4 can be written in matrix notation as 8, = Af

(5)

where 8, = e t ( x j ) , f = f (Vi),and A is a matrix with elements aj,i = O(Ui,Xj) AU. The ill-posed character of eq 2 is synonymouswith the fact that the matrix A is near singular and A should be inverted using a special technique in order to obtain a unique solution of A-l. Equation 5 is essentially the starting point for all three numerical methods. To calculate f, each of the methods uses a different mathematical technique to obtain A-l. Once a “stable” A-l is known f can be calculated. CAESAR. In CaesarWfi the problem of the near singularity of matrix A has been overcome by replacing A by a nonsingular matrix A,, using SVD. As a result eq 5 is replaced by 8, = A,f

(6)

in which A, is constructed by the decomposition

A, = UQ,Vt

(7)

U and V are orthogonal matrices (mXm and n x n , respectively) and Q, is a mXn diagonal matrix of the singular values of A,:

= diag(ql,qz,~3,...,~r,~l,...,~l) (8) The singular values qi are placed in decreasing order. Due to the fact that A is not exact only data up to qr should be trusted, where r is the numerical rank of A, defining the number of linearly independent simultaneous equations existing. Replacing values q+l, ...,qn by the largest singular value q1 (arbitrarily chosen large value) leads to Qp’ As a result A, is nonsingular. Equation 5 can now be solved with respect to f using a standard method. In order to obtain only positive values of f, use is made of Qp

which c is the iteration number and g j is themean deviation due to the error in e t ( l c j ) and e(U,Xj). Demanding the summation to be smaller than m would be asking for more accuracy than is inherently present in the data. In each iteration step r is increased by 1until the value is found satisfying eq 10. The computer program CAESAR offers three options for the local isotherm: the Hill-de Boer (HB) equation,’s a closed form of the virial equation,ls and the FrumkinFowler-Guggenheim (FFG) equation.12J3 In the testa described below CAESAR is used without modifications. PTHRegularization and REMEDI. Principal in the PTH regularization is the extra condition that f must be smooth. The PTH regularization minimizes the sum of squared second derivatives of f ( V ,on the energy interval A. Using this condition in combination with eq 5 leads to a new matrix equation:

(AtA + yH)f = A%, (11) where At is the transpose of matrix A, H is a constant nXn matrix representing the smoothness condition, and y is a non-negative Lagrange multiplier which determines the influence of matrix H on the solution. The term H serves as a nonsingular correction to the near-singular matrix AtA. Any conventional procedure to solve the set of linear equations may now be applied to find f, provided y is taken large enough. House46#47 used a Gaussian elimination algorithm with iterative improvement. Values of y were taken arbitrarilybut with the restriction that the obtained f(W was reasonably smooth. Disadvantages of the regularization as applied by House are that f ( V ,can become negative and that y is selected without a well-defined criterion. In order to overcome the disadvantages of the PTH regularization, we have modified the technique and introduced the NNLS proceduresg to find non-negative values for f only. A consequence of the use of NNLS is that smaller values of y are required to obtain an acceptable solution off than when a conventional procedure is used. For the selection of y a procedure similar to that for the selection of r in CAESAR is used. Starting with a lower limit of y, the vector f, is obtained and from that et,,. The solution is checked using eq 10. To calculate each new value of y , the bisection method is used. We will denote (59) Laweon, C. L.; Hanson, R. J. Solving Leust Squares Problems; Prentice-Hall: Englewood Cliffs,NJ, 1974.

2598 Langmuir, Vol. 9,No. 10, 1993

Koopal and Vos

U/RT U/RT U/RT Figure 1. Comparison of the calculated distribution functions with the true distribution function (solid curve) in the absence of an experimental error and in the presence of random errors in Ot of 1% and 5 % ,respectively. The three calculation methods and the magnitude of errors are indicated in the figure.

this modification of the PTH regularization as REMEDI, regularization method to obtain energy distributions. CAEDMON and CAEDMON-W. Although in CAEDMON the ill-posed-ness is not accounted for explicitly, it has been shownm that the requirement to obtain a solution free of spurious oscillations is that the number of knot points n on the energy distribution has to be chosen sufficiently lower than the number of points m of the experimental isotherm. As a result eq 5 becomes overdetermined and can be solved using a “least-squares” method. Sacher and Morrison& obtained the best results with the NNZS.59 The number of points of the energy distribution was selected arbitrarily, but n < m. Specific aspects of CAEDMON are that a discrete energy distribution is obtained and that only one local isotherm equation, viz., the closed form of the virial equation, can be used. In order to overcome the limitations of CAEDMON, we have made the following modifications: (1)f(U is seen as a continuous function, (2) the number of points, n (n < m),of the energy distribution function is selected on the basis of the errors in e&) and B(U,x),(3) the same types of local isotherm functions can be used as with CAESAR. In order to choose the right value of n, the solution is checked comparing with Bt as with CAESAR. The initial value of n is chosen equal to m, the number of data points. In each iteration step n is decreased by 1until the number (n1 2) is found satisfying eq 10. We will indicate this method as CAEDMON-W (Wageningen). (60)Papenhuijzen, J.; Koopal, L. K. In Adsorption from Solution; Ottewill, R. H., Rochester, C. H., Smith, A. L., Eds.;Academic Press: New York/London, 1983;p 211.

Table I. Mean Relative Deviation a(&) (%) as a Funotion of the Experimental Error in t+ for Each of the Three Methods Used A&(%) CAESAR REMEDI CAEDMON-W 0.0 0.002 0.002 0.003 0.1 0.03 0.04 0.02 1.0 0.28 0.28 0.26 5.0 1.39 1.48 1.31 10.0 3.41 3.67 2.78

Results and Discussion Test Isotherms. To illustrate the effectiveness of the methods outlined above, the inversion of some model gas adsorption data, generated by direct numerical evaluation of eq 1 for a HB equation as the local isotherm and a bi-Gaussian distribution function, is examined. The test function is described in more detail in ref 55. Apart from this bi-Gaussian distribution several other distribution functions have been briefly investigated. Results are comparable with those reported below. In Figure 1the original and the calculated distribution functions as obtainedwith the three methods are compared for the case that there is no error in the data or the local isotherm equation used. In all cases the result is very good. The absolute mean relative deviation, 6, between 4 and m

obtained for each method is shown in Table I. The results for the different methods are very similar. In practice experimental errors will always be present, and to simulate this situation Gaussian distributed random

Langmuir, Vol. 9, No. 10, 1993 2599

Adsorption on Heterogeneous Surfaces Table 11. Calculated Normalization Capacity as Percentage of the True Normalization Capacity AOt(%) CAESAR REMEDI CAEDMON-W 0.0 99.9 99.7 99.6 0.1 99.8 99.6 99.4 1.0 99.8 98.1 47.8 5.0 97.0 92.8 91.4 10.0 96.3 85.1 83.5 Table 111. Magnitude of the Characteristic Parameter Controlling the Instability of f(v) as a Function of the Experimental Error in Ot CAESAR AOt

0.0 0.1 1.0 5.0 10.0

r 21 13 8 8 4

REMEDI 'y

1.0 x 1.0 x 8.0X 2.0 x 2.5 X

1w lo-' lo-' 10-2 1 W

CAEDMON-W n

21 16 11 10 10

errors have been superimposed on the values of et. No error is added to the values of x j . In this way the standard deviation of the experimental isotherm was varied from 0.1 to 10%. The calculated distribution functions for a 1%and a 5% error are shown in Figure 1together with the true distribution. The observed values of 6(Ot) are summarized in Table I. The accuracy of the normalization capacity is shown in Table 11. When the adsorption is limited to monolayer coverage, as in the present case, the is directly related to the specific monolayer adsorption, rm, surface area. Hence, Armcorresponds to the error in the specific surface area based on rm.The results clearly indicate that the distributions obtained with CAESAR, REMEDI, and CAEDMON-W are very similar. CAESAR appears to give slightly better results for Armin the case of large experimental errors. From Table I it follows that 6 is increasing with an increasing experimental error. This is an inherent property of the present calculation methods. When large experimental errors are present in Bt eq 10 is satisfied, even when et and are not very close. A better fit would lead to details in f ( v) which are not corroborated by the experimental results. In Table I6 is noted as &et) in order to indicate that the deviation is due to an error in et. In Table I11the parameters r (CAESAR),y (REMEDI), and n (CAEDMON-W) are presented. The optimum values for r, y, and n change considerably by going from ABt = 0.0% to ABt = 0.1 % ,indicating that already a very small error in ABt leads to spurious oscillations in f (v). Table I11 also illustrates that the use of CAEDMON-W becomes impractical when the error is large because the number of points off(v) which can be calculated becomes very small. This is the most serious disadvantage of CAEDMON. An advantage of REMEDI relative to CAESAR is that y can vary smoothly, whereas r changes by 1 in each calculation cycle. As a consequence CAESAR may show small fluctuations on f ( v ) which are not observed with REMEDI. Compare for instance the region in between the two peaks in Figure 1. Another type of error, which in principle will always be present, is the implicit error in the chosen model for the local isotherm. In contrast to the experimental error of which mostly at least an a priori guess can be made, the error due to the choice of the local isotherm is difficult to quantify. To investigate this type of inaccuracy, a series of test calculations similar to those described in ref 55 were performed with local isotherm functions differing from the true one. For all methods the results are very similar to those

reported before.ws5 A typical result, calculated with REMEDI, is shown in Figure 2. Overestimatingthe lateral attraction (Figure 2a) results in a flattening of the curves and a small shift of the peak positions to lower energy values. This is due to the fact that the assumed local isotherm reaches higher coverages than the true local isotherm for a given value of x . When the lateral attraction is underestimated (Figure 2b), the peaks are sharper and shift to slightly higher energy values than the true peaks off( v). The reason for this behavior is that the assumed local isotherm increases less with x than the true local isotherm. When a local isotherm is used which has a symmetry that is different from the symmetry of the true local isotherm, extra peaks may appear in the distribution function. This is, for instance, the case by using the FFG equation while the correct local isotherm is a HB equation, compare Figure 2c. In Table IV the mean relative deviations 6 as observed under the various circumstances are summarized together with the relative deviation of the normalization constant, Arm. To indicate that the deviations are due to errors in the local isotherm, the parameters are denoted as 6(@ and Again it may be concluded that CAESAR, REMEDI, and CAEDMON-W give very similar results. From Table IV it follows that by increasing the lateral attraction the values of 6(e) decrease gradually. However, for the correct local isotherm function a minimum occurs in S(t9). In practice this minimum is very difficult to establish because of the experimental errors present and the approximate nature of the local isotherm. It should be realized that 6(e) is due to systematic errors. For a certain local isotherm these differences are minimized in the calculation, but they cannot disappear. In general, systematic deviations between the shapes of the calculated and the experimental overall isotherms will give an indication of the errors due to the wrong choice of For an adequate incorporation of this aspect in the numerical calculations further research is required. A wrong isotherm equation also leads to an error in rm. Table IV shows that Armcan be substantial when a wrong isotherm equation is used. Hence, comparison of rm calculated from the heterogeneity analysis with rmobtained in another way may provide a selection criterion for the local isotherm. Comparing 6(e) with (Tables I and IV) shows that for accurately measured isotherms 6(e) easily outweighs 6(Bt). In other words, for an accurate data set a large value of 6 is a strong indication that the chosen local isotherm is inadequate. Analysisof ExperimentalIsotherms. Nzon Carbon Blacks. To test the heterogeneity analysis with real adsorption data, surfaces of graphitizable carbon blacks are well suited. By heating such carbons in an inert atmosphere, the carbons graphitize, thereby creating a less heterogeneous surface, the result depending on the applied temperature. Graphitization at 2700 O C results for Sterling and Spheron carbon (Cabot Corp.) in homogeneous or nearly homogeneous surfaces.61 Annealing at lower temperatures will lead to partial graphitization. Accurately measured low-pressure N2 gas adsorption isotherms at 77.5 K are available for Spheron, Spheron 800 (heated for 12 h at 800 "0, Graphon (Spheron 27001, and Sterling MT 2700.62 These isotherm data have been investigated in several papers.41-45*56i57*62 There is evidence that the adsorption of N2 on Graphon and Sterling MT (61) Smith, W. R. Encyclopedia of Chemical Technology, 2nd ed.; Intsrecience: New York,1964,Vol. 4,p 243. Schaeffer,W. D.; Smith, W. R.; Polley, M. H. Znd. Eng. Chem. 1963,45, 1721.

2600 Langmuir, Vol.9,No.10,1993

Koopal and Vos I v) .

I

UiRT 0

I. Lo

cu

UIRT U/RT Figure 2. Effect of a wrong choice of the local isotherm in the adsorption energy distribution calculated with REMEDI. The true distributionfunction is indimted by the solid surve. The points represent the calculated distributions. (a) Overestimating the lateral attraction, the HB equation is used with an interaction parameter of 6.75 instead of 5.5. (b) Underestimatingthe lateral attraction, an interaction parameter of 3.9 is used instead of 5.5. (c) The FFG equation with w = -4.1 is used instead of the correct equation.

Table IV. Mean Relative Deviation a(@) and the Relative Deviation of the Normalization Constant AI’,(@) Wrong Choice of the Local Isotherm. -. .

..

HB. w’ = 5.5

HB; o’= 6.75 HB, w’ = 3.9

virial eq (Nz)

CAESAR

REMEDI

0.002 0.16 2.08 1.39 0.87

0.002 0.13 3.66 2.36 0.89

(both 46) Due to a

CAECDMON-W 0.003 0.11 2.08 1.34 0.89

FFG, w = -4.1 0 The HB (w’ = 5.5) equation is the correct local isotherm.

2700 is mobile.62fe3Hence, both the HB and the virial equation for Nz with coefficient(s) derived from 3D properties of Nz can be used as the local isotherm.s6 By using ideal 2D coefficients, no allowance is made for gassolid perturbations. For the analysis the surfaces will be regarded as patchwise heterogeneous. A correction for multilayer adsorption has not been applied to the adsorption data, but only adsorption values below the BET monolayer value have been used. Initial estimates of the energy interval used in the calculation are based on the condensation approximation. The final interval is chosen in such a way that the closure problems are small for all three methods. (62) Wewon, S.P. Adsorptiue Energies of CalorigenicPlesiomorpha of Graphite. Ph.D. Thesis, RensaelaerPolytechnic Institute, Troy, NY,

1975. (63) Van Dongen, R. H. Physical Adsorption on Ionic Solids. Ph.D. Theeis, Technical University Delft, The Netherlands, 1972.

CAESAR

REMEDI

CAEDMON-W

-0.1 -3.1 +13.4 +31.3 -21.7

-0.3 -4.8 +10.7 +33.4 -24.5

-0.4 -4.8 +10.5 +32.7 -24.6

First the total (average) error, which we have to consider in the calculations, is estimated. This error is composed of a part due to the choice of the local isotherm and a part resulting from the experimental error in Be. The error due to the inappropriateness of the local isotherm was obtained by applying CAESAR with a rank equal to the number of data points and the CA energy interval.& The affinity distribution functions for the four carbons, calculated with the three methods, using the virial equation for Nz as the local isotherm, are presented in Figure 3. Table V gives an overview of the parameters of the calculation. For each carbon the results obtained with the different methods are close: peak positions and magnitudes of the peaks are similar. The resolution obtained with CAESAR is slightly higher than that obtained with REMEDI or CAEDMON-W. The inherent problem of CAEDMON-Wis that for a smalldata set only a few points off( v) can be calculated. The mean relative

Adsorption on Heterogeneous Surfaces

0

9

CAEDMON-W

Langmuir, Vol. 9, No. 10, 1993 2601 REMEDl

CAESAR

f(U) 0

Y

0

9 0

U/RT U/RT UiRT Figure 3. Calculated adsorption energy distribution functions of NZadsorbed on a series of carbon blacks using the virial equation for N P as the local isotherm (patchwise heterogeneity): (a) Spheron; (b)Spheron 800; (c) Graphon (Spheron 2700); (d) Sterling MT 2700. The calculation methods me indicated at the top. Table V. Carbon Blackr-Na: Parameten Used in the Calculations of f l v ) Uring the Virial Equation er the Local Isotherm and the Resulting Valuer of b (Equation 12) and the Normalization Capacity rm P " 6 (% 1 rm (as (STP)/g) Spheron (m = 39, U h = 9.ORT; U,, = 21.0RT) CAEDMON-W n = 18 0.30 24.6 REMEDI y = 3 1 X le 0.32 26.0 CAESAR r=25 0.33 26.7 Spheron 800 (m = 60,U h = 9.5RT, U,, = 19.5RT) 0.52 25.2 CAEDMON-W n = 40 0.55 25.1 REMEDI y = 1.9 X le CAESAR r = 38 0.64 25.5 Graphon (m = 62, U h = 12.0RT, U- = 18.5R79 CAEDMON-W n = 12 1.15 29.0 REMEDI y = 0.7 X lo-' 1.21 29.9 CAESAR r = 20 1.18 30.2 Sterling MT 2700 (m = 49, U h = ll.ORT, U,, = 19.ORT) CAEDMON-W n = 20 0.90 3.29 1" y-16X le 1.36 3.26 0.96 3.29 CAESAR r = 27

deviation 6 between the measured and calculated isotherms (see eq 12) is in all cases very small. Both 6 and the normalization capacity, rm,hardly depend on the method used (see Table V). The resulta as obtained for Spheron and Spheron 800 show that the main peaks which can be recognized in both spectra remain the same. Heating to 2700 OC narrows the width of the distribution considerably in favor of the two peaks around 13RT. For SterlingMT 2 7 0 these two peaks

can no longer be distinguished. Note that the scale of the y axis for Figure 3a,b is different from that of Figure 3c,d. The agreement between the spectra for the different carbons and the similarity of the spectra as obtained with the different methods give faith to the calculated distributions and to the way they are calculated. In order to test the sensitivity of the methods for the choice of the local isotherm, the analysis is repeated with the HB (Nz)isotherm. Results are presented in Figure 4. As before the distributions obtained with CAEDMONW,REMEDI, and CAESAR are very similar for the same carbon. At the high-energy part of the distribution some instabilities can be present. Comparison of parts a-d of Figure 4 with their counterparts of Figure 3 shows that the distributions obtainedwith the HB (N2) equation tend to show more peaks than those calculated with the virial equation. This effect is very clear with Sterling MT 2700 where the HB equation results in a distribution function with two large peaks instead of one large peak as obtained with the virial equation. The two peaks are similar to thoee obtained for Graphon (using the HB equation), but somewhat narrower. The fact that the HB equation tends to give more detail than the virial equation is an indication that the virial equation is to be preferred as the local isotherm; compare Figure 2c. Also the fact that Sterling MT 2700 is known as a highly homogeneous surface is best illustrated with the distribution function obtained using the virial equation for N2 Comparison of older r e ~ u l t a ~ ~with ~ * $the ~ present ones

2602 Langmuir, Vol. 9, No. 10, 1993

Koopal and Vos

CAEDMON-W

n

I

REMEDl

CAESAR

U/RT

I

I

I

U/RT

U/RT

Figure 4. Adsorption energy distribution functions of Ns adsorbed on carbon blacks. The same as Figure 3 but with the HB equation for NZ (lateral attraction parameter 5.5) as the local isotherm.

shows that the presently derived distribution functions with the various methods are much closer to each other than those in previous comparisons. This is due to the fact that (1)the methods used are more similar (same continuous distribution function, non-negativity constraint, and local isotherm equation), (2) the applied boundary conditions (data range, energy interval) have been the same for the different methods, and (3) the same criterion is used (eq 10)tosuppressthe spurious oscillations in the distribution function. Chemically Modified Aerosils. Another set of data which has been used to calculatef (v) are the N2 adsorption isotherms of four chemically modified aerosils. The isotherms were measured gravimetrically at 78 K.@+j6 Preparation of the samples has been described by Berndt et alqMSample A is a heated Aerosil380, containing about one silanol group per square nanometer. Sample B is comparablewith sample A. The number of silanol groups is about the same as that of system A, but B contains imide groups instead of siloxane groups. In samples C and D the silanol groups as occurring in, respectively, systems A and B have been replaced by OSi(CH3)sgroups. The maximum adsorption values for the different samples range from 80 to 100%of the BET monolayer capacities. (64) H o w , W. A.; Jaroniec, M.; Bratier, P.;Fink,P. Thin Solid Film 1981,86,87. (66) House, W. A.; Jaroniec, M.; Bratier, P.; Fink,P. Thin Solid Film 1982,87,323. (66)Bemdt, H.; Born, G.;Braaer, P.;Fink, P.;Henneberg, K. H.;

Toparkus,H.Wbe. 2.Friedrich-Schiller-Uu. Jena, Math.-Natunuass. Reihe, 1981, 30, 1.

Just like House et al.@*@ we assume that the surface of the samples is patchwise heterogeneous and the adsorbate is mobile, so that both the virial and HB equatiolis for N2 as the local isotherm can be used to analyze the dataeS6 The data have been used without multilayer correction. The experimental error in the data has not been specified by House et al., but judging from the plotted isotherms the errors are small. The appropriateness of the local isotherm equations was checked by applying CAESAR and putting the rankequal to the number of data points, using the CA energy interval. The values for b(B) are small (about 0.06% of r,) for both the virial and the HB equation. Combination of this error estimate with an arbitrarily chosen experimental error in r of 0.25% provides the overalluncertainty for the calculation of the energy distribution. Results of the calculations using the virial equation for NBas the local isotherm are presented in parts a-d of Figure 5 for systems A-D, respectively. Each figure contains three distribution functions corresponding to the different solution methods. Table VI summarizes the parameters and the upper and lower limits of the energy intervals used in the calculations. Figure 5a shows the results for system A. The distribution functions obtained with CAESAR and FUMED1 are almost equivalent. For CAEDMON-W the number of points of the distribution function is much smaller, but the agreement with the distributions obtained with the other methods is fair. In Figure 6b f ( v ) for system B is

Adsorption on Heterogeneous Surfaces

Langmuir, Vol. 9, No. 10, 1993 2603 CAESAR

CAEDMON-W

UlRT

UlRT

UlRT

Figum 5. Calculated adsorption energy distribution functions of N2 adsorbed on a series of chemically modified aerosils. The virial equation for N2 is used BB the local isotherm (patchwise heterogeneity). The calculation methods are indicated at the top; panela a 4 represent the different silicas: (a) sample A; (b) sample B; (c) sample C; and (d) sample D. Table VI. Chemically Modified Aerosils-N,: Parameters

Used in the Calculations of 4V)Using the Virfal Equation as the Local Isotherm and the Resulting Values of I (Equation 12) and the Normalization Capacity rm Par6 (% 1 rm (mdg) A. SiOSiOH (m= 42, U,,,i,,= 6.0RT, U,, = 17.5RT) 0.24 67.6 CAEDMON-W n = 15 0.26 66.6 REMEDI y = 0.39 X 10-8 0.22 67.5 CAESAR r=14 B. SiNHSiOH (m= 40, U,,,i,,= 5.0RT, U,, - 18.0RT) n = 16 0.31 69.2 CAEDMON-W 0.33 72.0 REMEDI y = 0.20 X 10-9 0.33 67.7 CAESAR r = 19 C. SiOSI(CH& (m= 27, U- = 5.ORT, U- = 17.0RT) 0.08 90.9 CAEDMON-W rt = 12 0.13 116.7 REMEDI y = 0.20 x 10-8 0.05 90.0 CAESAR r = 15 D. SiNHSi(CH& (m= 24, U- = 5.ORT, U,, = 17.0RT) CAEDMON-W n 12 0.17 107.3 REMEDI y = 1.6 X 10-8 0.16 98.3 CAESAR r = 14 0.13 69.8

shown. Again the agreement between f ( V)obtained with CAESAR, REMEDI, and CAEDMON-W is very good. In Figure 5c,d the results for systems C and D are shown. Much less detail is present in the distribution functions calculatedwith REMEDI and CAEDMON-Wwith respect to both the result of CAESAR and the results shown for samples A and B. The smoothness may be due to the small number of experimental points, but it may also have

to do with the fact that the bulky trimethylsiloxanegroups screen the interactions with the rest of the surface. In order to show the effect of the local isotherm, the results obtained with the HB equation for Nz are presented in Figure 6. Again the results obtained with the different methods are similar. For system A a low affinity tail appears in the REMEDI spectrum, which does not occur with CAEDMON-W and CAESAR. System B shows relatively much detail, whereas systems C and D show little detail. In general a fair agreement exists between results obtained with the HB equation and with the virial equation. As with the carbon blacks, the HB equation leads to finer structure in the spectra. This is due to the fact that the smoothing is weaker with the HB equation than with the virial equation. In general, we may have underestimated the experimental error somewhat. The normalization capacities obtained by using the virial equation are considerably higher (about 80% ) than those for the HB equation. This may indicate that the lateral attraction is underestimated in the virial equation (see also Table IV). Both rm(V)and I',(HB) are larger than the BET value of I',.a*=This indicates that a multilayer correction should be applied to the data. The rising part of f(U)at low values of UIRTalso points in this direction. The results shown in Figures 5 and 6 compare fairly well with those obtained by House et al.a@ using HILDA and CAEDMON (NNLS). However, especially the HILDA results show finer structure which is probably unrealistic.

2604 Langmuir, Vol. 9, No. 10, 1993

Koopal and Vos

CAEDMON-W

CAESAR

t

?

U/RT

U/RT

U/RT

Figure 6. The same as Figure 5, but the HB equation for Nz is used as the local isotherm.

Concluding Remarks The LIA methods, used in combinationwith generalized cross validation plus physical constraints to preprocess the data, are well developed to calculate an approximate distribution function. The error in the data determines whether the CA (AS1) or the LOGAl (AS2, RJ) method should be used. When the error in the data is of the order of a few percent, the CA method has to be used. The LIA methods are especially suited for the heterogeneity analysis when only limited information is known (window of data, neglect of lateral interactions, poorly known local isotherm, etc.). The calculated distribution function may, for example, be used to decide whether or not one of the existing analytical overall isotherm equation@ is suited to describe the data.37 The comparison of the results of the numerical methods shows that regularization in combination with practical constraints for the selection of the stability or smoothing parameter (as used in CAESAR and REMEDI) is very well suited to calculate the adsorption energy distribution function. Although also CAEDMON-W may give good results, it has the disadvantage that the method is less rigorous and that only a limited number of points of a spectrum can be calculated. The distribution functions obtained are not absolute. The amount of detail in a distribution depends on both the quality of the data and the quality of the local isotherm. (67)Van R i e d i j k , W.H.;Koopal, L. K. InEnuironmental Particles; Buffle, J., van Leeuwen, H.P., Eds.;Lewis Publishers: London, 1992; Vol. 1, Chapter 12.

Spectra obtained for the graphitizable carbonsshow a clear correspondence, and the results obtained hardly depend on the analysis technique for a given local isotherm equation. The same applies to the spectra obtained for the different silicas. These results indicate that the obtained distribution functions are free of large spurious peaks; most details shown are a direct image of significant inflections of the isotherm. Both CAESAR and REMEDI are well suited as tools to study heterogeneity on a routine basis. They could, for example, be of help in production or quality control of activated carbonsor other particle surfaces. By comparing not only BET areas and porosities, but also the N2 affinity spectra of different batches of a certain product, valuable additional information is attained. Affinity spectra allow a detailed comparison of different samples, even when the peaks cannot be related directly to chemical information. In the last decades the emphasis has been on the development of methods which deal with the ill-posedness of the problem. With the “regularization” methods (PTH, REMEDI, CAESAR, INTEG)fil55*” and the FT technique5*this problem is well solved. The challenge for the future is to improve the selection of the “smoothing” parameter. One way to do this for the regularization methods is to take advantage of the achievements made with the LIA methods.3m7 Applying a smoothing spline routine in combination with GCV plus the physical constraint of a positive first derivative of e&) to the isotherm data allows us to obtain a fairly good estimate of the error in the data. The possible error due to the use of a particular local

Adsorption on Heterogeneous Surfaces

Langmuir, Vol. 9, No. 10,1993 2606

isotherm can subsequently be investigated by considering method.69 In this case the number of points of the distribution function (nknot points) should be sufficiently the spline result as the accurate isotherm and applying lower than the number of experimental points, m, in order CAESAR (REMEDI) with a rank (smoothing parameter) to obtain a stable solution. The best possible solution for determined by the machine precision. The standard a specified experimental error is obtained by using eq 10. deviation of the spline combined with the estimated error The initial value of n is chosen equal to m, and in each due to the use of a particular local isotherm can now be iteration step n is decreased by 1until eq 10 is satisfied. used to find the final rank (smoothing parameter) in an In this way peaks due to experimental error are smoothed objectiveway, and the affinity spectrum can be calculated. away. To investigate whether the obtained distribution A personal computer would be suited for such calculations. function is stable, one can still further decrease n by 1or A thus-obtained result could be compared with the LIA 2 and observe how strongly this affects the summation of result which relies entirely on presmoothing of the data. eq 10and the number of peaks in the distribution function. The constraint presented by the use of a specific local So far the discussion has been limited to the methods isotherm in REMEDI or CAESAR is replaced in the LIA used in the adsorptionliterature. A challenge for the future method by further constraintsto the data points depending is to screenthe physical scienceliterature for other methods on the actual LIA. for inverting noisy linear operator equations. Especially For Langmuir- or FFG-type local isotherm equations it the CONTIN program of P r o v e n ~ h e deserves r ~ ~ ~ ~ atten~ is also possible to apply the smoothing spline GVC routine tion. Similarly as with CAESAR and REMEDI the with three additional constraints and to use the thusCONTIN program allows for specific constraints and knowledge about the experimental errors to “stabilize” obtained spline as an accurate representation of the data the solution. The advantage of Provencher’s method is from which the distribution function can be calculated that it is a mature method in which the underlying directly with REMEDI or CAESAR without further mathematics have been worked out more satisfactorily iterations (smoothing parameter or rank determined by machine precision). Similarlyas with the LIA m e t h 0 d s ~ ~ 3 ~ than with the methods presented here. Comparison of the results of Provencher’s CONTIN method applied to this procedure allows for the calculation of confidence the data sets used here with the present results will intervals for the distribution function. This method is certainly be of help for the further development of methods relatively quick, and its result can be compared with that to establish the optimal distribution function for a given obtained in the “regular” way. data set. The presented methods work with three options for the local isotherm. This restriction is not essential: when (68) Olivier, J. P.; Conklin, W. B. Presented at the Effects of Surface Heterogeneity in Adsorption and Catalysis on Solids Symposium. another local isotherm equation is preferred, the chosen (69) Lastoskie,C.; Gubbm, K. E.; Quirke, N. Langmuir 1993,9,2693. isotherm equation can, in principal, be incorporated in (70) Provencher, S. W. Comput. Phys. Commun. 1982,27,213. any of the regularization methods. If the local isotherm (71) Provencher, S. W. Comput. Phye. Commun. 1982,27,229. (72) Although the acronym originally refers to the computer program is fairly steep and only known numerically (see e.g. refs used for the calculation off( 0, we will use it also to indicate the method 68 and 69)) the simplest solution is to solve the integral as a whole. The same applies to the other methods. equation for the distribution function by using the NNLS (73) One of the reviewers made ua aware of the existence of CONTIN.