Density Functional Theory Analysis of the Influence of Pore Wall

Lourdes F. Vega. 2007 ... Testing the feasibility of using the density functional theory route for pore size distribution calculations of ordered micr...
12 downloads 0 Views 211KB Size
Langmuir 2002, 18, 6845-6856

6845

Density Functional Theory Analysis of the Influence of Pore Wall Heterogeneity on Adsorption in Carbons Suresh K. Bhatia* Department of Chemical Engineering, The University of Queensland, Brisbane QLD 4072, Australia Received February 21, 2002. In Final Form: June 17, 2002

Density functional theory for adsorption in carbons is adapted here to incorporate a random distribution of pore wall thickness in the solid, and it is shown that the mean pore wall thickness is intimately related to the pore size distribution characteristics. For typical carbons the pore walls are estimated to comprise only about two graphene layers, and application of the modified density functional theory approach shows that the commonly used assumption of infinitely thick walls can severely affect the results for adsorption in small pores under both supercritical and subcritical conditions. Under supercritical conditions the Henry’s law coefficient is overpredicted by as much as a factor of 2, while under subcritical conditions pore wall heterogeneity appears to modify transitions in small pores into a sequence of smaller ones corresponding to pores with different wall thicknesses. The results suggest the need to improve current pore size distrubution analysis methods to allow for pore wall heterogeneity. The density functional theory is further extended here to allow for interpore adsorbate interactions, and it appears that these interaction are negligible for small molecules such as nitrogen but significant for more strongly interacting heavier molecules such as butane, for which the traditional independent pore model may not be adequate.

Introduction Carbon-based microporous materials, such as activated carbons, carbon molecular sieves, and carbon fibers, as well as mesoporous materials, such as carbon xerogels, have wide applications as adsorbents in separations and as supports in catalysis. Apart from their traditional use in water purification, there has also been a long-standing interest in the use of carbons in protective gas masks amid filters for a variety of applications, such as in mines, warfare, and accidental radioactive gas release. More recently there has been considerable growth in their applications in gas-phase separations for industrial pollution abatement, in air separation, and in other industrial processes. In addition, newer applications such as adsorptive storage of natural gas in carbons for transportation applications continue to emerge.1,2 Consequently, there is considerable interest in the understanding and analysis of adsorption in such materials, with a large and growing body of literature on the subject. The analysis and mathematical models of the adsorption equilibrium also provide important tools for characterization of the carbon, as the adsorption equilibrium is strongly affected by adsorbent structure and properties.3 Among the key properties of the adsorbent affecting the adsorption equilibrium is its heterogeneity, by virtue of which there are differences in adsorption energy between different sites or regions of the solid surface. In the modeling of adsorption equilibrium, this is generally considered in terms of the distribution of adsorption energies,4,5 or a distribution of pore sizes with the fluid* Fax: +61-7-3365 4199. Telephone: +61-7-3365 4263. E-mail: [email protected]. (1) Quinn, D. F.; MacDonald, J. A. Carbon 1992, 30, 1097. (2) Matranga, K. R.; Myers, A. L.; Glandt, E. D. Chem. Eng. Sci. 1992, 47, 1569. (3) Gregg, S. G.; Sing, K. S. W. Adsorption, Surface Area and Porosity; Academic Press: New York, 1982. (4) Ross, S.; Olivier, J. P. On Physical Adsorption; Interscience: New York, 1964.

solid interaction or adsorption energy being a function of pore size.5,6 The latter formulation leads to the well-known generalized adsorption isotherm (GAI)

Fj(P) )

∫0∞F(P,H) f(H) dH

(1)

in which f(H) is a pore volume distribution function for pores of size H, P is the bulk gas-phase pressure, and F(P,H) is the local isotherm in a pore of size H. For carbons the slit-pore model is commonly used, with the pore size H interpreted as the corresponding slit width. The slit model serves as a useful basis for characterization, though it is clearly approximate, as it is unlikely that all graphitic domains in carbons are parallel. Micropores within the domains are, however, expected to be parallel, given the layered structure of the domains. For the local isotherm there exists a variety of models in the literature based on traditional analyses that consider the adsorbate in a pore as a single uniform phase, thereby overlooking density gradients, as well as newer approaches that consider the molecular texture and nonuniformity of the adsorbed phase. Among the former class of models are those7-9 based on the well-known Dubinin-Radushkevich and DubininAstakhov models, those10 based on the Langmuir isotherm, and a more recent vacancy solution theory formulation.11-13 Among the molecular modeling based local isotherms are those using density functional theory6,14-16 or molecular simulation techniques.17,18 An important feature inherent to the majority of the local isotherm models referred to above is the assumption (5) Rudzinski, W.; Everett, D. H. Adsoprtion of Gases on Heterogeneous Surfaces; Academic Press: New York, 1992. (6) Lastoski, C.; Gubbins, K. E.; Quirke, N. Langmuir 1993, 9, 2693. (7) McEnaney, B. Carbon 1988, 26, 267. (8) Bhatia, S. K.; Shethna, H. K. Langmuir 1994, 10, 3230. (9) Ismadji, S.; Bhatia, S. K. Langmuir 2001, 17, 1488. (10) Jagiełło, J.; Schwarz, J. A. Langmuir 1993, 9, 2513. (11) Bhatia, S. K.; Ding, L. P. AIChE J. 2001, 47, 2136. (12) Ding, L. P.; Bhatia, S. K. Carbon 2001, 39, 2215. (13) Ding, L. P.; Bhatia, S. K. AIChE J., accepted.

10.1021/la0201927 CCC: $22.00 © 2002 American Chemical Society Published on Web 08/10/2002

6846

Langmuir, Vol. 18, No. 18, 2002

Bhatia

of infinitely thick pore walls and the associated Steele19 10-4-3 potential when the adsorptive is LennardJonesian (LJ). This potential has the form

[(

) ( )

2 σsf φsf(z) ) A 5 z

10

σsf z

4

-

σsf4

]

3∆(0.61∆ + z)3

(2)

where A ) 2πFssfσsf2∆, z is the distance between an adsorbate molecule and the solid surface, ∆ is the interplanar distance in the carbon, Fs is the LJ site density of the solid, and sf and σsf are fluid-solid LJ parameters following the Lorentz-Berthelot rules. The above expression for the potential at a distance z from the surface approximates that due to infinite stacked planes separated by the spacing ∆. In reality, however, most carbons have pore walls with only a few, typically one to three, planes of carbon atoms in a hexagonal close packed structure. Under this circumstance it may be anticipated that the potential in eq 2 overestimates the strength of the fluid-solid interaction and therefore overpredicts the amount adsorbed. However, the extent of the error introduced has never been examined. The quantification and investigation of the error are particularly important because the 10-4-3 potential in eq 2 is routinely used in density functional theory (DFT) or simulation based models,14-17 as well as other more conventional approaches,12,13,20 in estimating pore size distributions from adsorption data. For example, nitrogen or argon adsorption isotherms are commonly used to determine pore size distributions in conjunction with an appropriate isotherm model; however, the error in the capillary coexistence curve due to the assumption of infinitely thick pore walls, and the consequent overestimation of pore size are unknown. A further complicating feature is that the carbon framework enveloping the pores is highly disordered,21 as a result of which the wall thickness may be expected to be random and nonuniform. In a recent attempt to improve on the infinite wall thickness assumption, Ravikovitch et al.22 have taken one of the two walls of a pore to comprise a single layer, while keeping the other wall infinitely thick. However, while such an arbitrary assignment of wall thickness may be adequate in special cases, there is a need for a more general formulation that allows for the structural disorder and randomness in the wall thickness in a less confining way. This article reports on a modified density functional theory formulation that permits a random distribution of wall thicknesses to be considered. It is also shown that the mean wall thickness can be readily evaluated from the pore volume distribution and that nonuniformity in the wall thickness can be considered without invoking any new parameters or constants. The density functional theory procedure used here is the modification23 of the Denton-Ashcroft (DA) procedure,24 which is computa(14) Olivier, J. P.; Conklin, W. B.; Szombathely, M. V. Stud. Surf. Sci. Catal. 1994, 87, 81. (15) Seaton, N. A.; Walton, J. P. R. B.; Quirke, N. Carbon 1989, 27, 853. (16) Ravikovitch, P. I.; Vishnyakov, A.; Russo, R.; Neimark, A. Langmuir 2000, 16, 2311. (17) Davies, G. M.; Seaton, N. A. Langmuir 1999, 15, 8736. (18) Vishnyakov, A.; Ravikovitch, P. I.; Neimark, A. V. Langmuir 1999, 15, 8736. (19) Steele, W. A. Surf. Sci. 1973, 36, 317. (20) Ismadji, S.; Bhatia, S. K. J. Colloid Interface Sci. 2001, 244, 319. (21) Stoeckli, H. F. Carbon 1990, 28, 1. (22) Ravikovitch, P. I. Carbon’01, International Conference on Carbon, Lexington, July 14-19, 2001. (23) Bhatia, S. K. Langmuir 1998, 14, 6231.

tionally simpler than others, and has been validated23 against simulation results. The extent of the error in the isotherms due to assumption of infinitely thick pore walls is also investigated for both supercritical and subcritical conditions, and the impact on pore volume distributions thereby determined is discussed. Although we utilize the modified23 Denton-Ashcroft approach here, it is to be noted that several alternate density functional theory-based models exist. Insightful reviews of the various methods have been provided by Evans25 and Davis,26 and we shall not detail these here. However, we do note that the most popular of these is that due to Tarazona27 and is commonly used in determining pore size distributions6,14-16 or modeling single component adsorption.28 We have chosen the modified DA procedure because of its computational advantage and because it is also applicable to the binary adsorption case. On the other hand, the Tarazona prescription has yet to be extended to binary systems. For the latter case the Kierlik-Rosinberg29 procedure has found some application30,31 but is computationally more intensive than the modified DA method of the author.23 Nevertheless, it should be noted that the choice of DFT procedure does not affect the generality of the formulation, which relates to the consideration of pore wall nonuniformity and could be applied to any other prescription to obtain comparable results. Model Development Density functional theory is adapted here to permit a distribution of pore wall thicknesses. Since the DFT approach is now well established, with several variants reviewed in the literature,25,26 we focus our attention on the DA method24 used here and also employed elsewhere by the author.23 An important contribution is the determination of the mean pore wall thickness from pore structure information, and its incorporation in the theory. In what follows we first develop the relation between the pore size distribution and the wall thickness in slit-pore carbons, before illustrating its incorporation in the DFT approach. Pore Wall Thickness in Carbons. Key to the present development is the concept of a random pore wall thickness, as may be expected for the highly disordered carbons that are considered appropriate as models for activated carbons.30 Although the infinite wall thickness assumption is common to almost all adsorption models, as indicated above, it is readily estimated that the actual pore wall thickness is likely to be very small. To this end we consider a porous carbon with a micropore surface area Sg. Recognizing that this area corresponds to that of the two surfaces bounding the pore walls, we may write

Sg )

2N0 γFsMc

(3)

where N0 is the Avogadro number, γ is the mean number of graphene layers per wall, Fs is the number of carbon (24) Denton, A. R.; Ashcroft, N. W. Phys. Rev. A 1991, 44, 8242. (25) Evans, R. In Fundamental of Inhomogeneous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992. (26) Davis, H. T. Statistical Mechanics of Phases, Interfaces and Thin Films; Wiley-VCH: New York, 1996. (27) Tarazona, P. Phys. Rev. A 1985, 31, 2672. (28) Tan, Z.; Marconi, V. M. B.; Gubbins, K. E. J. Chem. Phys. 1989, 90, 3704. (29) Kierlik, E.; Rosinberg, M. L. Phys. Rev. A 1991, 44, 5025. (30) Rodriguez-Reinoso, F.; Linares-Solano, A. Chem. Phys. Carbon 1988, 21, 1. (31) Lebowitz, J. L. Phys. Rev. 1964, 133, A895.

DFT Analysis of Adsorption in Carbons

Langmuir, Vol. 18, No. 18, 2002 6847

atoms per unit area in a plane, and Mc ()12) is the atomic weight of carbon. Equation 3 provides

γ)

2N0 FsMcSg

(4)

which conveniently expresses the mean number of layers in the walls in terms of surface area. The surface site density for the graphene planes of the crystallites of which carbons are comprised, Fs, has the value 0.3817 Å-2, so that, on substitution of the values of N0 and Mc, eq 4 yields

γ)

2629.9 Sg

(5)

in which Sg has the units square meter per gram. Most microporous carbons used as adsorbents typically have surface areas in the range 800-1500 m2/g, though higher values are also reported in the literature. It is now readily estimated that for most carbons the mean number of layers lies in the range 1.75-3.25, which suggests that existing adsorption models severely overestimate the adsorbatecarbon interaction by virtue of the assumption of infinitely thick pore walls. Further, eq 5 sets the upper limit of the surface area of slit-pore carbons with an underlying crystallite structure as 2629.9 m2/g, as the pore walls must comprise at least one graphene plane. The mean number of graphene layers in the pore walls can also be related to the pore size distribution by virtue of the relation

Sg )

∫0∞H2 f(H) dH

(6)

where f(H) is the pore size distribution, defined such that f(H) dH represents the pore volume per unit mass in micropores of width lying in [H, H + dH]. Equations 4 and 5 yield

γ)

h N0H FsMcvµ

(7)

the nonuniform pore size distribution is one manifestation of the heterogeneity, it is to be expected that there is also considerable randomness in the pore wall thickness. The latter is of much significance in view of the small mean wall thickness, as discussed above. Since the fluid-solid potential in general has a range of a few molecular diameters, any random variations in wall thickness from two layers to one, or to three or four layers, can affect the potential field in the pore significantly. To allow for such possibilities, one may consider a probability distribution for the number of graphene layers in the pore walls. In general, such a wall thickness distribution would depend on the method of preparation of the carbon, but if we accept a model of disorder with no bias or central tendency in the wall thickness, the Poisson distribution for the number of layers arises as the natural choice. Accordingly, we may write

p(n) )

(γ - 1)n-1e-(γ - 1) , n ) 1, 2, ... (n - 1)!

where p(n) is the probability that the wall has n layers. The above model considers the fact that a wall must have at least one layer of carbon atoms, so that the number of additional layers (n - 1) is considered as a random variable with mean (γ - 1). Besides being the natural distribution for the case of complete or maximum disorder, the Poisson distribution is also the simplest, as it has only one parameter (the mean number of layers). Since this mean is related to the pore size distribution, it is evident that the present model does not introduce any new parameters in the point size distribution (PSD) estimation. While other distributions could be chosen, they have at least one additional parameter such as a standard deviation. Density Functional Theory. The density functional theory employed is that due to Denton and Ashcroft,24 subsequently modified23 for improved results in binary systems. Here we adapt the single component formulation to incorporate the structural disorder in terms of wall thickness distribution. For this we define the grand potential Ω[Flm] in a pore of width H, with the left wall having l layers and right wall having m layers, as

where vµ is the total specific micropore volume

vµ )

∫0∞f(H) dH

(8)

and

H h )

vµ ∞f(H)

∫0

H

(9) dH

is the harmonic mean pore width. Upon substituting the values of Fs, Mc, and N0 in eq 7, we obtain

0.1315H h γ) ) vµ

0.1315 ∞f(H) dH 0 H



(10)

in which vµ has the units cubic centimeters per gram, H h is expressed in angstroms, and f(H) is in cubic centimeters per gram angstrom. Equations 8-10 express the mean number of layers per pore wall in terms of the pore size distribution. Often the latter is the primary quantity obtained from structural characterization, with the surface area Sg being a derived quantity. Microporous carbons are inherently disordered materials,30 with considerable structural heterogeneity. While

(11)

Ω[Flm] ) F[Flm] +

∫dr Flm(r)[Vlm(r) - µ]

(12)

in which Flm(r) is the density profile while Vlm(r) is the external potential in the pore and µ is the chemical potential of the adsorptive. Further, F[Flm] is the intrinsic Helmholtz free energy, given by

F[Flm] ) Fid[Flm] + Fex hs[Flm] + 1 dr Flm(r) dr′ Flm(r′) φatt(|r-r′|) (13) 2





where Fid[Flm] is the ideal gas free energy at the local density Flm, following



Fid[Flm] ) kT dr Flm(r){ln[Flm(r)λ3] - 1}

(14)

and Fex hs[Flm] is the hard sphere excess with the last term on the right-hand side of eq 13 representing the attractive contribution in the fluid-fluid potential. The parameter λ represents the de Broghie wavelength of the adsorptive at temperature T, and k represents the Boltzmann constant. Equation 13 ignores interactions among adsorbate molecules in neighboring pores, that must arise because the wall separating these is no longer assumed to be infinitely thick. However, it will be subsequently

6848

Langmuir, Vol. 18, No. 18, 2002

Bhatia

shown, on the basis of a mean field approximation, that such interactions can in general be neglected for disordered carbons. Unconstrained minimization of the grand potential in eq 13 provides the Euler-Lagrange relation

Flm(r) ) Fb exp[-βVlm(r) + c(1)(r;[Flm]) - c0(1)(Fb) -





att

in which β ) 1/kT, and we have used the condition that Flm(r) f Fb for the bulk fluid for which Vlm(r) ) 0. Here c(1)(r;[Flm]) is the one-particle direct correlation function (DCF) of a hard sphere system, given by

δFex hs[Flm] c(1)(r;[Flm]) ) -β δFlm(r)

c(1)(r;[Flm]) ) c0(1)(Fjlm(r))

(17)

where

∫dr′ F(r′) w(|r-r′|;Fjlm(r))

(18)

The normalized weight function w(r;F0) is obtained from

w(r;F0) )

∂c(1) 0 /∂F0

(19)

which ensures that the first functional derivative of the one-particle DCF in eq 17 with respect to the density yields the correct two-particle DCF c(2)(r;F0) in the uniform limit (i.e. F(r) f F0). In the present work, as in the original DA treatment24 and the earlier study of the author,23 the wellknown Perens-Yevick approximation result31 for c(2)(r;F0) is used. The result for the weight function in the onedimensional slit geometry has been provided in the earlier works.23,24 Potential Models. A Lennard-Jones (LJ) 12-6 pair potential

φ(r) ) 4

12

6

[(σr) - (σr) ]

) -; r < rm

(21) (22)

of the LJ potential, where rm ()21/6σ) is the location of the minimum of this potential. For the adsorbate-carbon interactions the LJ parameters are obtained using the Lorentz-Berthelot rules. where ff and cc are the well-

fc ) xffcc

(23)

σfc ) (σff + σcc)/2

(24)

depths for the the fluid-fluid adsorptive and carbon-

(25)

where u(r) is the repulsive part of the LJ potential. In evaluating the potential field in carbon slit pores it has in the past been common to use the Steele 10-4-3 potential,19 based on infinitely thick pore walls. In the present case the potential field due to a wall with n layers is obtained by summing the 10-4 potential due to each layer,19 to obtain n-1

φfs(z,n) ) 2πFsσfs2fs

[ ( ) ( )] σfs

2

∑ i)0 5 z + i∆

10

-

σfs

z + i∆

4

(26)

where z is the distance between an adsorbate molecule and the solid surface, and ∆ is the interplanar spacing. For carbon the value of ∆ is 3.4 Å, and cc/k has the value 28 K.19 The total potential in a pore of width H is now obtained by adding the contributions from each wall, based on eq 26, to obtain

Vlm(z) ) Φfs(z,l) + Φfs(H-z,m)

(27)

Heterogeneity. Structural heterogeneity is considered in the present study through the presence of a pore size distribution as well as a wall thickness probability distribution. The density profile obtained upon solution of the single pore equations presented above can now be averaged over the pore cross section as well as the wall thickness distribution to obtain the mean adsorbate density in a pore of width H as ∞

Fγ(H) )

∑ m)1



1

p(l) ∫0 Flm(z) dz ∑ H l)1

p(m)

H

(28)

where the profile is assumed to be one-dimensional and the thicknesses of the opposing pore walls are assumed to be uncorrelated. When a pore size distribution is also considered, eq 28 can be averaged over the pore volume to obtain the mean absorbate density in the pore structure as

(20)

with well depth  and size parameter σ is used for the intermolecular interactions in this work. The attractive part of the pair potential is obtained according to the Weeks-Chandler-Anderson division

φatt(r) ) φ(r), r > rm

∫0∞[1 - exp(-u(r))] dr

(16)

In the DA24 formulation of the DFT, to account for nonlocal effects, the one-particle DCF for the inhomogeneous fluid is approximated by that of a homogeneous fluid at a suitably weighted density Fjlm(r), that is,

c(2)(r;F0)

d)

att

β dr′ Flm(r′) φ (|r-r′|) + βFb dr′ φ (r′)] (15)

Fjlm(r) )

carbon potentials, respectively, and σff and σcc are the corresponding size parameters. The adsorptive hard sphere diameter for evaluating the DCFs is estimated by means of the Barker-Henderson expression

〈Fγ〉 )

∫0∞Fγ(H) f(H) dH

1 vµ

(29)

where 〈Fγ〉 represents the mean density. For the calculations reported here a gamma distribution of micropore volume is considered, following

q(b+1)Hbe-qH f(H) ) vµ Γ(b+1)

(30)

in which the parameters b and q are related to the mean pore size, H0, and relative standard deviation of the PSD, s, following

b)

(1 - s2)

q)

s2 1 s H0 2

(31) (32)

DFT Analysis of Adsorption in Carbons

Langmuir, Vol. 18, No. 18, 2002 6849

Figure 1. Variation of mean number of layers in a pore wall with relative standard deviation in the PSD, for total pore volume 0.5 cm3/g.

where

H0 )

∫0∞Hf(H) dH

1 vµ

(33)

and

s)

[

]

∫0∞(H - H0)2f(H) dH

1 1 H0 vµ

1/2

(34)

Thus, the three parameters H0, s, and vµ completely specify the PSD. In terms of these parameters the mean number of layers in the pore wall is readily obtained following eqs 10 and 30-32, as

0.1315H0(1 - s2) γ) vµ

Figure 2. Variation of total amount adsorbed with dimensionless bulk density, for various values of mean pore width. Potential parameters correspond to those for methane. The primary figure displays the isotherm in the low-pressure range, while the inset shows the isotherm over a larger pressure range. Solid lines represent the isotherms for the disordered carbon with heterogeneous pore walls, while the dashed lines represent the isotherms for infinitely thick walls.

or more layers essentially offered the potential of an infinite wall. Consequently, the summations in eq 28 were taken only up to five terms, with the modification 4

(35)

Figure 1 depicts the variation of the mean number of layers in the pore walls, γ, with relative standard deviation s, for vµ ) 0.5 cm3/g and various values of the mean pore size H0. Most carbons have mean micropore width ∼8-10 Å and a value of s ∼ 0.2. As seen in Figure 1, for such carbons the value of γ is in the vicinity of about 2. It is also evident that wider pore size distributions lead to thinner pore walls, which is also easily deduced from eq 35. Clearly, any pore wall must comprise at least one graphene layer, and as is seen in Figure 1, this limit is attained for s lying in the range 0.7-0.8 for typical carbons with H0 ≈ 8-10 Å. This provides an upper limit for the dispersion in the PSD using the gamma function form. Results and Discussion The above model has been applied in the present work to study the effect of pore wall heterogeneity on single component equilibria under supercritical as well as subcritical conditions. The carbon micropore is idealized as an infinite slit of width H (defined as the center-tocenter distance between the carbon atoms on the opposing pore walls), so that planar symmetry may be assumed. Consequently, the density and other spatial variables are only a function of the normal coordinate z, measured from one wall. The integral equations following from eq 15 were then discretized in a uniform grid of mesh size σff/10, which yields sufficient accuracy. The discretized equations were solved iteratively, mixing old and new solutions to promote convergence. Further, in estimating the mean adsorbate density for a carbon with mean number of layers in the walls, γ, following eq 28, it was found that walls with five

p(5) ) 1 -

∑ p(n)

(36)

n)1

Adsorption under Supercritical Conditions. For supercritical conditions the LJ parameters of the adsorptive were taken to be those of methane, for which σff ) 3.81 and ff/k ) 148.2 K. Initially, the effect of temperature on the adsorption isotherm was studied, for a carbon having a gamma distribution of PSD with micropore volume 0.5 cm3/g, mean pore width H0 ) 2.25σff, and relative standard deviation in PSD s ) 0.25. For such a carbon, eq 35 provides the mean number of layers in the pore walls as γ ) 2.25. Under supercritical conditions the dimensionless temperature T* () kT/ff) must have a minimum value of about 1.29, following the well-known approximation ff/kTc = 0.78 for LJ fluids.32 Figure 2 depicts the computed dimensionless isotherms for methane in the carbon with the above PSD, at various dimensionless temperatures. In this figure the solid lines represent the isotherms for the above disordered carbon with finite pore walls having γ ) 2.25, while the dashed curves provide the corresponding results assuming instead infinitely thick pore walls and the Steele 10-4-3 potential, as is the usual practice. The figure shows the isotherms in the law density, or Henry’s law, region while the inset depicts the isotherms in the region approaching saturation. As seen in the figure, while the difference between the two isotherms is small (within 5%) at high bulk densities, there is a considerable deviation (≈20-100%) in the Henry’s law region. This deviation is particularly large at lower temperatures (i.e. T* f 1.33). Clearly, at low pressures the adsorption in (32) McQuarrie, D. A. Statistical Mechanics; University Science: Sausalito, CA, 2000.

6850

Langmuir, Vol. 18, No. 18, 2002

Bhatia

Figure 4. Effect of relative standard deviation in the PSD on the relative overall isotherm for a carbon with mean pore width 2.25σff and total pore volume 0.5 cm3/g, at T* ) 1.75.

Figure 3. Variation of relative amount adsorbed, Finf/Fγ, with dimensionless bulk density at dimensionless temperature T* ) 1.75 (a) in a pore of width 2.25σff at various values of mean number of layers in pore walls, γ, and (b) in pores of various sizes at fixed mean number of layers in pore walls γ ) 2.25.

real carbons, with finite pore walls, is considerably weaker than that anticipated by the infinite wall assumption. This is because in the low-pressure region the adsorbateadsorbate interactions are weak, and the fluid pore wall interaction dominates. At high pressures, however, the interaction between adsorbate molecules is very significant and the effect of modification of the fluid-wall interaction potential is much less severe. Figure 3a depicts the effect of mean pore wall thickness, presented by the parameter γ, on the variation of the ratio Finf/Fγ with dimensionless bulk density in a pore of width 2.25σff at T* ) 1.75. Here Finf is the adsorbate density in the pore, assuming infinitely thick walls, with the Steele potential expression in eq 2 being therefore used. The ratio approaches a constant in the Henry’s law region, as is to be expected. In this region the ratio is approximately 2 for γ ) 1.5, indicating that the assumption of infinite pore walls leads to an overestimation of the Henry’s law

constant for the adsorption by this factor. At γ ) 3 this overestimation drops to a factor of about 1.25. With increase in bulk density the error initially increases slightly, before dropping rapidly at a dimensionless bulk density Fσff3 of about 10-4. This initial increase in error is due to cooperative effects, whereby the larger adsorption in the infinite wall case leads to enhanced adsorbate interactions that serve to attract more adsorbate molecules. Beyond a dimensionless bulk density of about 10-4 (or bulk pressure of about 0.065 bar), the adsorbateadsorbate interactions become important even for the finite wall case and the ratio rapidly drops to near unity as the effect of wall interaction decreases in relative significance. Figure 3b depicts the effect of pore size on the variation of the ratio Finf/Fγ with dimensionless bulk density Fσff3 at T* ) 1.75 and γ ) 2.25. The overestimation in the Henry’s law constant by the infinite wall assumption increases with decreasing pore size, being about 60% higher for H ) 2σff. Even for H ) 3σff (or pore width 11.43 Å) use of the infinite wall assumption leads to Henry’s law constant being higher by about 35%. The low-pressure region is one of much importance because of the interest in removal of trace organics and pollutants. Consequently, the present results would suggest that reliable process design in this region must consider the disorder of the carbon and its effect on the adsorption. Figure 3b also shows that the overestimation in the low-pressure region shows a much more marked and abrupt increase for the pore width 3σff, at a dimensionless bulk density of about 10-3, compared to other pore sizes. At this pore size it is evident that two layers of adsorbate can be accommodated, while the smaller pore sizes 2σff and 2.5σff can accommodate only one layer, and the larger size 4σff can accommodate three layers. As has been shown in an earlier article,23 the packed density under supercritical conditions is the highest in pores of sizes of about 3σff, where exactly two layers can be accommodated. Consequently, adsorbate-adsorbate interactions are stronger and cooperative effects most significant at this pore width, leading to the more prominent peak in the ratio Finf/Fγ. Figure 4 depicts the effect of the relative standard deviation in the PSD on the overall isotherm, for a carbon with mean pore width 2.25σff and micropore volume 0.5

DFT Analysis of Adsorption in Carbons

cm3/g, at T* ) 1.75. The value of γ in each case is calculated on the basis of eq 35 and lies in the range 1.69-2.22 for the values of s used in Figure 4. The pattern is similar to the results in Figure 3 for individual pores, with the ratio 〈Finf〉/〈Fγ〉 approaching a constant in the Henry’s law region and falling to a near unity value at high pressures. However, it is evident that the overestimation of the adsorbate density and the Henry’s law constant increased with an increase in the dispersion of the PSD. This is to be expected on the basis of eq 35, which predicts a decrease in mean wall thickness with increase in relative standard deviation s. Even for the value s ) 0.25 (for which γ ) 2.11) it is seen that the Henry’s law constant is overestimated by over 60% when the infinite wall assumption is used. All of the above results strongly indicate the importance of considering realistic pore wall thicknesses for carbons in adsorption modeling. While we have investigated the effect of finite pore wall thickness on supercritical adsorption above, it is to be noted that the PSD, on which the mean pore wall thickness depends [cf. eq 10], is itself usually determined from low-temperature adsorption isotherms using adsorptives such as nitrogen and argon. The present results suggest that existing procedures for PSD determination using regularization or other methods10,14-17 must be modified to consider finite wall thickness. To demonstrate the need for this, we examine below the effect of finite wall thickness on low-temperature adsorption isotherms of nitrogen. Adsorption under Subcritical Conditions. Lowtemperature adsorption using adsorptives such as nitrogen or argon at their boiling points is perhaps the most widely used method for pore structure characterization. In this method eq 1 is inverted, using the measured isotherms along with DFT or other models to provide the local isotherm, to determine the PSD. The present adaptation of the DFT to incorporate finite mean pore wall thickness with a random distribution in number of layers provides a new local isotherm, and here we examine this feature considering nitrogen at 77 K as the adsorptive in carbon micropores. The actual implementation to the determination of the PSD is, however, consigned to a future of study. To apply the present DFT to nitrogen adsorption, initially the model was first used to extract suitable LJ parameter values by fitting data on the variation of bulk vapor pressure as well as saturated liquid- and gas-phase densities with temperature. This yielded σff ) 3.45 Å and ff/k ) 102.626 K for nitrogen, which is only marginally different from the values of σff ) 3.572 Å and ff/k ) 93.98 K obtained by Lastoskie et al.6 using the Verlet and Weiss34 approximation for the effective hard sphere diameter. Indeed, when the latter approximation was used in place of eq 25, the parameter values obtained were closer to those reported by Lastoksie et al. However, here we used the above values based on use of the more exact BarkerHenderson expression in eq 25. These parameters were then used in the DFT to predict adsorption isotherms in pores of different sizes. To obtain the isotherms, the adsorption and desorption scanning curves for any pore size and wall thickness combination were obtained by solving the model equations at gradually increasing and decreasing pressures, respectively, starting with an empty pore in the former case and a nearly saturated pore in the latter case. This provides the two minima in the grand potential, with a phase transition when the grand (33) Bhatia, S. K. Chem. Eng. Sci. 1998, 53, 3239. (34) Verlet, L.; Weiss, I. J. Phys. Rev. A 1972, 5, 939.

Langmuir, Vol. 18, No. 18, 2002 6851

Figure 5. Predicted nitrogen isotherms in pores of various sizes, for carbon with infinitely thick pore walls. Solid lines represent present calculations, while symbols represent the molecular simulation results of Lastoskie et al.6

potentials at these two solutions are equal. From these solutions the isotherm was determined from the points of lowest grand potential, at any pore size and wall thickness combination, and the wall thickness averaging performed following eq 28. Figure 5 depicts the predicted nitrogen isotherms at 77 K (solid lines) in pores of various sizes, as well as the molecular simulation results of Lastoskie et al.,6 for the infinite wall case, showing good correspondence. Here H* is a dimensionless pore size defined as H/σff and P0 is the vapor pressure. To facilitate the comparison, the potential model of eqs 26 and 27 was used with n ) 5, since the Steele 10-4-3 potential in eq 2 was employed for the Monte Carlo simulations of Lastoskie et al. In addition, the carbon LJ parameters used were the same as those of Lastoskie et al., having the values σcc ) 3.416 Å and cc/k ) 30.14 K. These are only slightly different from those used in the computations discussed above. For estimating the pressure at any density, the following bulk equation of state, obtained from the generic free energy functional, was used.

Fb2 P ) kTFb + 2

∫bulkφatt(r) dr - kTFbC(1) 0 (Fb) + F kT∫0 C(1) 0 (F) dF b

(37)

The DA prescription for the DFT used here matches the molecular simulation results very well, though for very low pore sizes of 2.5σff and smaller a small deviation appears, as shown in Figure 5. While the Tarazona27 formulation used by Lastoskie et al.6 gave slightly better results at the lowest pore size, it is evident that the present approach provides an adequate basis for studying the effect of disorder in the wall thickness on the isotherms. Figure 6 depicts the adsorption isotherms of nitrogen at 77 K in carbon pores of various sizes, in the presence of disorder with γ ) 2.0. The dashed curves for all pore sizes correspond to the isotherms when both pore walls are infinitely thick, while the doted curves represent the isotherms when both walls have only one graphene layer. These present the limiting isotherms at any pore size, and in the presence of pore wall disorder, the isotherm for

6852

Langmuir, Vol. 18, No. 18, 2002

Bhatia

is, however, somewhat less, as seen in Figure 6b for H* ) 5 and H* ) 12, and at large pore size it has virtually no effect on the transition. This is because in such pores surface tension effects dominate over the wall-fluid interaction, and the well-known Kelvin equation limit is approached. Effect of Pore-Pore Interactions. The above results have been based on an independent pore model, in which the adsorbate density profile obtained from eq 15 for any combination of wall thicknesses does not depend on the adsorption in neighboring pores. However, when the walls separating these are thin (for example having only one graphene layer), there is concern about the significance of interactions between adsorbate molecules in such pores. The impact of such interactions may be expected to be most significant under subcritical conditions, whence phase transitions occur. Consequently, we have examined this effect here for the low-temperature nitrogen isotherms discussed above. The effect of pore-pore interactions is most conveniently incorporated in the present formulation by modifying the external potential in eq 27 to include those in the model for the grand potential in eq 12. Thus, the external potential in a pore of size H with the left wall having l layers and the right wall having m layers is now expressed as

Vlm(r) ) Vext lm (r) + 1



∑ p(l′)∫0 [∫Fl′l(r′,Hl) φ(|r-r′|) dr′]h(Hl) dHl + 2l′)1

1





∑ p(m′)∫0 [∫Fmm′(r′,Hm) φ(|r-r′|) dr′]h(Hm) dHm 2m′)1 ∞

(38)

Figure 6. Comparison between predicted nitrogen isotherms for carbons having pore wall disorder with γ ) 2.0 (solid lines) and those for carbons having infinitely thick walls (dashed lines). Dotted lines represent the isotherms for pores with both walls having only one layer: (a) isotherms for H* ) 2.5 and (b) isotherms for H* ) 3, 5, and 12.

any pore wall thickness combination, as well as the effective isotherm obtained from eq 28, must lie between these limits. Thus, the space between these curves represents the domain of possible isotherms at any pore size in the disordered carbons. The solid curves in Figure 6 correspond to the effective isotherms, based on eq 28, for γ ) 2.0. From the figure it is evident that disorder has a profound effect on the isotherms, and the effective isotherm no longer displays the distinct monolayer transitions in the narrow micropores (H* ) 2.5 and 3). Instead the isotherm appears almost continuous with several smaller transitions corresponding to the pores with various wall thicknesses. From Figure 6 it is evident that except for larger pore sizes and high relative pressures the infinitely thick pore wall assumption can severely overestimate the amount adsorbed. Consequently, current methods for PSD determination may be expected to somewhat overestimate the pore size in the carbon. As in the case of supercritical adsorption, considerable deviation of the isotherm is evident at low and moderate pressures, which is the important region for micropore size determination. The effect of disorder on the capillary condensation transition

where Vext lm (r) represents the fluid-solid interaction component in eq 27, while recognizing the planar symmetry in the infinite slit pore model assumed. Here we assume the existence of a mean density field Flm(r,H) in a pore of size H with left and right walls having l and m layers, respectively, in the ensemble of pores in the microporous structure. The first summation on the righthand side of eq 38 represents the interactions with the pore to the left having random size Hl with probability distribution h(Hl) given by

h(Hl) )

f(Hl)/Hl

∫0 [f(H)/H] dH ∞

(39)

Further, the left pore shares a common wall with the target pore of size H, which becomes the right wall of the former pore. However, its left wall has a random number of layers l′ with probability distribution p(l′) following eq 11, so that the mean field density in this pore is Fl′l(r,Hl), and the summation over l′ averages over this distribution. In a similar way, the last term on the right-hand side of eq 38 represents the interactions with the pore on the righthand side of the target pore of size H. Since the consideration of pore-pore interactions requires that we abandon the independent pore concept, we now define the overall grand potential in the structure as ∞

ΩT[F] )



p(l) ∑ p(m)∫0 Ω[Flm(r,H)]h(H) dH ∑ l)1 m)1 ∞

(40)

which averages the grand potential in eq 12 with respect to the pore size as well as left and right wall thicknesses.

DFT Analysis of Adsorption in Carbons

Langmuir, Vol. 18, No. 18, 2002 6853

The grand potential ΩT[F] may now be minimized with respect to the density field Flm(r,H) to obtain an EulerLagrange equation similar to eq 15:

Flm(r,H) ) (1) (1) Fb exp[-βVext lm (r) + C (r;[Flm]) - C0 (Fb) -





β dr′ Flm(r′,H) Φatt(|r-r′|) + βFb dr′ Φatt(r′) ∞

β

∑ F(l′)∫0 ∫Fl′l(r′,Hl) Φ(|r-r′|) dr′ h(Hl) dHl l′)1



β



∑ p(m′)∫0 ∫Fmm′(r′,Hm) Φ(|r-r′|) dr′ h(Hm) dHm] m′)1 ∞

(41) As expected, the density profile in a pore of size H is now correlated with the profile in pores of other sizes and wall thicknesses in the structure, and it depends on the pore size distribution. The above formulation of the model including interactions considers only the effect of adjacent pores, which is sufficient, as fluid-fluid interactions have influence only over a few molecular diameters. In a recent study, Patrikiejew et al.35 have also considered such interactions, but rather than using our mean field approach for the structure, they examine particular configurations of wall thicknesses with the periodicity of the lattice. The present formulation considers all possible configurations with wall thickness probabilities following eq 11, and it is therefore well-suited for determining the adsorption in the entire structure. To study the effect of pore-pore interactions, eq 41 has been solved here for nitrogen adsorption at 77 K using the carbon and nitrogen LJ parameters given above and the iterative technique indicated earlier. For simplicity, we have chosen a structure of uniform pore size, but with a wall thickness distribution following eq 11 with γ ) 1.75. As before, the adsorption and desorption scanning curves were obtained by solving the model at gradually increasing and decreasing pressures, respectively, starting with an empty structure in the former case and one nearly saturated with adsorbate in the latter case. In the presence of interpore interactions, however, it is necessary to conduct a detailed study of the phase behavior to unequivocally determine the stable solutions and the isotherm, as opposed to the independent pore model where the scanning curves define the two stable (or metastable) minima of the grand potential. Indeed, in a recent lattice gas study of adsorption in highly disordered porous solids Kierlik et al.36 report the scanning curves to represent a sequence of metastable states without an underlying equilibrium transition between them. Instead, a complex free energy landscape is observed with the adsorption and desorption branches enveloping a large number solutions, from among which an equilibrium curve could be constructed by joining the points of minimum grand potential. Thus, because of the influence of interpore interaction, the isotherm is not simply determined by selecting the points of lower grand potential at any pressure from the two solutions for any pore wall combination determined by the present scanning procedure. Consequently, we first compare here the scanning curves with and without the consideration of interpore interactions. (35) Patrikiejew, A.; Reszko-Zygmunt, J.; Rzysko, W.; Sokołowski, S. J. Colloid Interface Sci. 2000, 228, 135. (36) Kierlik, E.; Monson, P. A.; Rosinberg, M. L.; Sarkisov, L.; Tarjus, G. Phys. Rev. Lett. 2001, 87, 055701.

Figure 7. Effect of interpore interaction on nitrogen adsorption and desorption scanning curves, in the presence of pore wall disorder with γ ) 1.75, for structures having uniform pore size H* ) 2.5 and 3. Solid and dashed lines represent the adsorption and desorption scanning curves, respectively, without interaction, while dotted and dash-dot-dot curves represent the corresponding results considering interpore interaction: (a) overall scanning curves and (b) scanning curves for pores with both walls having only one layer.

Figure 7 depicts the scanning curves for nitrogen at 77 K in structures having uniform pore sizes H* ) 2.5 and H* ) 3. Part a depicts the overall scanning curves for the disordered structure, showing the effect of pore-pore interactions to be small, though for the larger pore size there is a significant effect on the adsorption curve at high relative pressures of about 10-4. In these figures the solid and dashed lines represent the desorption and adsorption curves, respectively, for the independent pore model without pore-pore interactions, while the dotted and dash-dot-dot lines represent the corresponding results in the presence of interactions. An interesting feature is the appearance of some hysteresis in the latter case for H* ) 2.5, for which there is no hysteresis in the absence of interactions (the solid and dashed curves coincide for H* ) 2.5). This is predominantly due to the effect of interactions on the fluid in pores having single layer walls, with pores having both walls with two or more layers behaving almost as in the independent pore model. This is evident from Figure 7b, which depicts the hysteresis

6854

Langmuir, Vol. 18, No. 18, 2002

behavior for pores with both walls being single layer, showing significant differences over the same relative pressure range as that for which differences appear in Figure 7a. In particular, we see the appearance of hysteresis effects for H* ) 2.5 in the presence of interactions, though they are absent in the independent pore case in pores with both walls being one layer thick. As another example of the effect of interpore adsorbate interactions, the model was applied to the case of more strongly interacting adsorptive molecules. For this the fluid LJ parameters were taken to be those corresponding to butane. Although the LJ model is only a crude approximation for a polyatonic molecule such as butane, it suffices to illustrate the trends with regard to the interpore fluid interactions. The LJ parameters for butane were estimated by simultaneously fitting its saturation vapor pressure as well as saturated gas- and liquid-phase densities at its normal boiling point of 272.4 K, which provided σff ) 4.8978 Å and ff/k ) 371.8677 K. These parameters as well as the carbon LJ parameters given above were then used in the solution of eq 41. As for the case of nitrogen adsorption, a structure having uniform pore size was assumed, but with a wall thickness distribution following eq 11 with γ ) 1.75. Figure 8 depicts the scanning results for butane at 272.4 K in structures of pore sizes H* ) 2 and H* ) 3. As for the case of nitrogen adsorption, a significant effect of interpore interactions is seen at H* ) 2 for the case of pores having only one layer in each wall (Figure 8b). However, even for the overall isotherm the effect of the pore-pore interactions is not negligible. In the relative pressure range 10-5 to 10-4, as much as 25-30% error occurs, predominantly because of the large interaction effect in pores with at least one single layer wall. In the case of H* ) 2, however, no hysteresis is present in both cases (with and without interactions). A large error (2530%) also occurs for H* ) 3 over certain pressure ranges, as the transitions occur at slightly lower pressures when interactions are considered. This pore size of about 14.7 Å is well above the mean micropore width of most activated carbons (typically 8-10 Å), and the results therefore suggest that pore-pore correlation effects are important for adsorption of large molecules such as butane or other bulkier species in carbons. To some extent the magnitude of this effect depends on the assumed pore wall thickness distribution. However, with other distributions that are peaked and have a greater central tendency, the effect will still be significant at low mean wall thickness. Given the similarity of the scanning curves and the weak (though not insignificant) effect of interactions, it seems likely that, as in the independent pore case, the equilibrium isotherms in the present case approximately correspond to the points of lower grand potential for each pore wall thickness combination. Thus, we anticipate that such points approximately correspond to the stable solutions, and other solutions between the scanning curves correspond to metastable states. While this procedure is exact for the isotherm in the case of independent pores, in the presence of interactions, phase transitions in pores of a given wall thickness combination must have some influence on the density in other pores. Nevertheless, this influence may be considered small because of the mild effect of interactions on the scanning curves, as seen in Figures 7 and 8, and approximate isotherms constructed ignoring this influence. However, a more comprehensive study of the solutions and phase behavior in the presence of interactions needs to be carried out to confirm this. Such a study is currently being planned in our laboratory and will be reported in due course. Nevertheless, at this

Bhatia

Figure 8. Effect of interpore interaction on butane adsorption and desorption scanning curves, in the presence of pore wall disorder with γ ) 1.75, for structures having uniform pore size H* ) 2 and 3. Solid and dashed lines represent the adsorption and desorption scanning curves, respectively, without interaction, while dotted and dash-dot-dot curves represent the corresponding results considering interpore interaction: (a) overall scanning curves and (b) scanning curves for pores with both walls having only one layer.

time it is instructive to compare the equilibrium isotherms in the absence of interactions with the curves determined as above for the interacting system. Parts a and b of Figure 9 depict the results for nitrogen for the structures of pore size H* ) 2.5 and H* ) 3, respectively. In this figure the dashed curves correspond to the isotherms for pores with both walls being one layer thick when the pore-pore interaction is neglected, while the dotted curves represent the corresponding result when the interaction is considered. Also shown in these figures are the overall isotherms, averaged over all combinations of wall thicknesses, with the solid lines representing the result with interactions neglected and the dash-dot-dot lines giving the result considering interactions. From the results it is seen that at the pore size 2.5σff pore-pore interactions are of significance for pores with both walls being only one layer thick. However, the overall isotherm is only marginally affected by the interactions. At the slightly larger pore size of H ) 3σff, it is seen that even

DFT Analysis of Adsorption in Carbons

Figure 9. Effect of interpore interaction on nitrogen isotherms, in the presence of pore wall disorder with γ ) 1.75, for structures having uniform pore size (a) H* ) 2.5 and (b) H* ) 3. Solid and dashed lines represent the overall isotherm and that in pores with both walls having only one layer, respectively, without considering interaction, while dotted and dash-dot-dot curves represent the corresponding results considering interpore interaction.

for pores with both walls being one layer thick the effect of interactions is small, though there is a small (≈5%) increase in density due to interactions after the monolayer transition. As in the case of the smaller pore size, the overall isotherm is only marginally affected with a very small increase in density at high pressures. From these results it would appear that for nitrogen the effect of porepore interactions is significant only at small wall thicknesses (single layer walls) and at small pore size, consistent with the observation of Patrikiejew et al..35 In actual carbons with a randomly distributed wall thickness, the resulting effect of the interactions on the overall isotherm appears to be small and may be neglected in comparison to the effect of the disorder itself. Figure 10 depicts the results obtained as above for butane, for the structures of uniform pore size H* ) 2 and H* ) 3. In this case, even at the larger pore size of H* ) 3 there is a significant effect of interaction on the isotherms for single walled pores, with the transition occurring earlier when interaction is considered. From the overall isotherms it is evident that large differences (as much as 40-50%) over narrow pressure ranges can occur, con-

Langmuir, Vol. 18, No. 18, 2002 6855

Figure 10. Effect of interpore interaction on butane isotherms, in the presence of pore wall disorder with γ ) 1.75, for structures having uniform pore size (a) H* ) 2 and (b) H* ) 3. Solid and dashed lines represent the overall isotherm and that in pores with both walls having only one layer, respectively, without considering interaction, while dotted and dash-dot-dot curves represent the corresponding results considering interpore interaction.

firming the effect of interactions to be more important for bulkier molecules such as butane, as discussed earlier. Conclusions In this research we have examined the effect of randomly distributed pore wall thickness on adsorption in activated carbons. For typical carbons it is seen that the mean wall thickness corresponds to only about two graphene layers, while current adsorption models generally assume infinitely thick walls. Density functional theory is adapted here for carbons with a random pore wall thickness, and it is seen that the infinitely thick wall assumption significantly overpredicts the adsorption under supercritical conditions, particularly in the Henry’s law region where errors as large as 100% occur. This error increases with decrease in mean wall thickness and with increase in dispersion of the pore size distribution. Under subcritical conditions it is seen that the isotherms in small pores are most significantly affected, with multiple small transitions

6856

Langmuir, Vol. 18, No. 18, 2002

arising from pores of the same size but having different wall thicknesses. The results suggest that current methods of pore size distribution analysis require further development to incorporate pore wall heterogeneity. The effect of interaction between adsorbate molecules in adjacent pores is also considered here, and the DFT has been extended to incorporate this feature in a mean field sense. It is seen that while the effect of interaction is small for nitrogen adsorption, it is significant for bulkier molecules such as

Bhatia

butane that interact more strongly with each other. These results demonstrate the weakness of the independent pore model for adsorption of large molecules. Acknowledgment. Financial support of this research by the Australian Research Council through the Discovery Grants scheme is gratefully acknowledged. LA0201927