Liquid-state theory of hydrocarbon-water systems: application to

J . Phys. Chem. 1992, 96, 8582-8594

8582

Llquld-State Theory of Hydrocarbon-Water Systems: Application to Methane, Ethane, and Propane Leo Lue and Daniel Blankschtein* Department of Chemical Engineering and Center for Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 021 39 (Received: May 7, 1992)

We have studied the structural and bulk thermodynamic properties of hydrocarbon (methane, ethane, and propane)-water systems as well as pure water using the site-site Omstein-Zemike (SSOZ) equation under a variety of different closure relations in order to compare the quantitative predictive Capabilities of the various closures. For the hydrocarbon-water systems, the simple point-charge (SPC) potential was used to model water, and the optimized potentials for liquid simulation (OPLS) were used to model the hydrocarbons. For pure water, predictions were also made for other water potential models. We solved the SSOZ equation with the hypemetted-chain(HNC) closure to determine the pair correlation functionsof water. We then analyzed the structural and bulk thermodynamic properties of methane, ethane, and propane at infinite dilution in water using various closure relations for the hydrocarbonwater pair correlation functions. We find that the HNC closure, which is the closure that has been utilized almost exclusively to predict bulk thermodynamic properties of interaction-site fluids, performs rather poorly. Specifically, we find that the HNC closure consistently underpredicts the magnitudes of both the solute partial molar volume and the solute-solvent interaction energy, grossly overpredicts the magnitude of the residual chemical potential, and gives the incorrect sign of the enthalpy of solution. On the other hand, we find that two recently developed closures, the Martynov-Sarkisov (Mas) and Ballone-Pastore-Galli-Gazzillo ( B E G ) closures, which have not been utilized so far in conjunction with the SSOZ equation, yield reasonable predictions of the structural and bulk thermodynamic properties of the hydrocarbon-water systems studied. In particular, utilizing the SSOZ-BEG equation, the predicted temperature variation of the residual chemical potential over the relatively broad range 5-80 OC was found to be in very good agreement with the experimental data. Note that the residual chemical potential is directly related to the Henry’s law constant, which, in tum, can be utilized to predict solubilities. In addition, we have developed an analytical expression for the residual chemical potential, appropriate for interaction-site fluids, in terms of pair and direct correlation functions at full coupling for the various closures examined in this paper. To date, an expression of this type was available only for the HNC closure. This new expression facilitates the calculation of the residual chemical potential by eliminating the previous need to perform a numerical integration over the coupling constant, thus making the computationof the chemical potential simpler and more efficient. Finally, we have also tested the accuracy of the equivalent-site approximation (ESA), a perturbation method which was developed by Curro and Schweizer to study long polymeric chains by treating all the sites in a given molecule as equivalent, on the predictions of the structural and bulk thermodynamic properties of propane at infinite dilution in water. Note that of all the n-alkanes, propane poses the most severe challenge to the ESA. Interestingly, we find that, already for propane, the ESA yields predictions of bulk thermodynamic properties which are within 5% of those obtained using the rigorous calculations.

1. Introduction The quantitative prediction of structural and bulk thermodynamic properties of hydrocarbon-water systems represents a challenging problem in liquid-state theory. These systems are particularly interesting because they are among the simplest to exhibit the hydrophobic effect.’ The central goal of this paper is to examine the ability of various closure relations, when combined with the site-site Ornstein-Zernike (SSOZ) equation? to provide accurate quantitative predictions of structural and bulk thermodynamic properties of the hydrocarbons methane, ethane, and propane in water from knowledge of their interaction potentials only. It is noteworthy that previous attempts to predict the bulk thermodynamic properties of interaction-site fluids of the type considered in this paper have utilized almost exclusively the hypemetted-chain (HNC) closure. As will be shown in this paper, the quantitative predictive capabilities of the SSOZ-HNC equation are rather poor. In view of this, it appears particularly important, both from a fundamental perspective and an applied perspective, to study whether the use of other closure relations can improve the quantitative predictive capabilities of the SSOZ equation beyond those of the HNC closure. For this purpose, we have solved the SSOZ equation with two recently developed closure relations: (i) the Martynov-Sarkisov (MaS) closure3 and (ii) the Ballone-Pastore-Galli-Gazzillo (BPGG) ~ l o s u r e .We ~ find that the SSOZ-MaS and SSOZBPGG equations perform significantly better than the SSOZHNC equation when applig to the hydrocarbon-water systems considered in this paper. In particular, we will show that utilizing *To whom correspondence should be addressed.

the SSOZ-BEG equation, the predicted temperature variation of the residual chemical potential over the relatively broad range 5-80 OC is in very good agreement with the experimental data. This result is significant, because,among other things, the residual chemical potential is directly related to the Henry’s law constant, which, in turn, can be utilized to predict solubilities. In addition to presenting a comprehensive comparative study of the predictive capabilities of various closure relations, we have also developed an analytical expression for the residual chemical potential, appropriate for interaction-site fluids, in terms of pair and direct correlation functions at full coupling for the various closure relations examined in this paper. To date, an expression of this type was only available for the HNC closure. The new expression facilitates the calculation of the residual chemical potential by eliminating the previous need to perform a numerical integration over the coupling constant. We hope that the availability of this new analytical expression will encourage the use of the recently developed closures in future studies of the thermodynamic behavior of molecular fluids. In the calculations proposed above, the number of equations that need to be solved increases approximately as the square of the number of distinct interaction sites belonging to the molecules in the system. Note that two interaction sites are considered to be distinct if the correlation functions associated with each one of them are different. In order to reduce the number of equations to a tractable number for molecules possessing many distinct sites, such as polymer chains, Curro and Schweizer developedS a perturbation method where, at zeroth order, all the sites in a given molecule are treated as equivalent. Note that two sites are considered to be equivalent if the correlation functions associated with each one of them are equal. In this paper, we refer to this

0022-365419212096-8582$03.00/0 0 1992 American Chemical Society

Liquid-State Theory of Hydrocarbon-Water Systems zerothsrder approximation as the equivalent-site approximation (EBA). Clearly, the ESA is exact in the case of ethane. For other n-alkanes, one expects the ESA to become better as the hydrocarbon chain length increases due to the decreasing importance of the terminal CH3 groups as compared to that of the internal CHI groups. In this paper, we will test the ability of the ESA to predict structural and bulk thermodynamic properties of propane at infinite dilution in water. Note that of all the n-alkanes, propane poses the most severe challenge to the S A . We will show that, interestingly, already for propane, the ESA yields predictions of bulk thermodynamic properties which are within 5% of those obtained using the rigorous calculations. Before we go on to describe the theoretical formulation, a brief summary of previous works which are particularly relevant to this paper is in order. The SSOZ equation has been used in the past to study some aspects of hydrocarbon-water systems. For example, by utilizing the water-water pair correlation functions obtained from X-ray scattering, in the context of the SSOZ equation with the Perm-Yevick (PY) closure, h t t and Chandler were able to make6 accurate predictions of the hydrophobic interaction free energy, as well as other thermodynamic properties, of various n-alkanes in water as a function of temperature. It is also noteworthy that the SSOZ equation with the HNC closure has been utilized to study pure water, as well as water in the presence of several types of solute species. Specifically, Pettitt and Rossky applied’ the SSOZ-HNC equation to study the structure of pure water using various potential models of water. Yu et al. utilized**9this theory to predict the solvation properties of various ions in water. Ichiye and Chandler usedl0 the SSOZ-HNC equation to study the propane-water system. The remainder of the paper is organized as follows. In section 2, we describe the central elements of the theory which is utilized in this paper. This includes (i) a description of the potential models of water and the hydrocarbons, (ii) a brief overview of the SSOZ equation and the various closures that we utilize, (iii) the development of a new analytical expression to calculate the residual chemical potential, and (iv) a few computational details. In section 3, we present the results of our predictions of the structural and bulk thermodynamic properties of hydrocarbons (methane, ethane, and propane) at infinite dilution in water and compare them with available experimental data and computer simulation results. In addition, we present the results of these predictions for propane in water using the ESA in order to assess the accuracy of this approximation. In section 4, we present a discussion of the results reported in section 3. In section 5 , we summarize the central findings of this paper. In addition, for completeness and for the benefit of those interested in quantitative predictions, in Appendix A we present new results of a detailed study of the structural and bulk thermodynamicproperties of pure water for various potential models of water. Finally, Appendices B and C contain the details of the derivations of some of the formulas utilized in the main body of the paper. 2. Theory

A. Potential Functim. We utilize the interaction-site model” to describe the interactions between molecules. In these models, sites are located at various positions in a molecule, typically on atoms or groups of atoms. The interaction energy between two molecules is taken to be the sum of the interaction energies between the sites in each molecule. In the interaction-site models that we utilize, the potential functions between sites are given as the sum of a long-ranged (1/r) Coulombic term and a short-ranged (12-6) Lennard-Jones-type potential; that is,

The Journal of Physical Chemistry, Vol. 96,No. 21, I992 8583 TABLE I: Intenetion Panmetera of Hydrocprbonsa z,M

0

CHI CH3 CHI

0 0

~___

Ca MI kcal f6/mol

8529 6935 5935

3167 2396 1674

See ref 12.

nardJones potential. The parametes Z,M, A , M ~ M and , C,M,Mare available in the literature for various atoms and groups of atoms (that is, for various aM). To estimate the parameters AaMdMt and CaMatM,,the following combining rules are utilized:I2 AaMa‘M‘ = (AaMaMAa‘M‘a‘M‘)

(2.2)

and (2.3) In the interaction-site models that we employ, the bond lengths and bond angles are held fixed. Consequently, the only manner in which a molecule can alter its conformationis through a rotation around its dihedral angles. Under these constraints, the molecules which we consider, namely, water, methane, ethane, and propane, are rigid. In the calculations of the structural and bulk thermodynamic properties of hydrocarbons in water, we have used the simple point-~hargel~ (SPC) potential model of water. Details of this potential model are given in Appendix A (see Table VII). In addition, for completeness and for the benefit of those interested in quantitative predictions, we have also calculated the structural and bulk thermodynamic properties of pure water using other potential models of water. Details of these calculations are also presented in Appendix A. To model the hydrocarbons, we have utilized the optimized potentials for liquid simulation (OPLS) developed by Jorgensen et a1.I2 The interaction parameters z, A, and C associated with this model for CH4, CH3, and CH2 are listed in Table I. For the OPLS,the carbon-carbon bond length is 1.53 A, and the CH3-CH2-CH3bond angle is 112O. B. Sitesite Omstein-Zernike Equation and Closure RehtioDn The Ornstein-Zernike (OZ) equation has had a great deal of success in the field of simple fluidsi4 Note that a simple fluid is defmed as a fluid consisting of molecules which interact through spherically symmetric potentials. Chandler and Andersen extended2J5the OZ equation to describe molecules which contain several interaction sites, that is, to interaction-site fluids. The resulting equation is called the sitesite Omstein-Zernike (SSOZ) equation and is written below in matrix form in momentum (k) space h(k) = G(k)t(k)G(k)+ @(k)t(k)ph(k)

(2.4)

where h is the total correlation function, C is the direct correlation function, G describes the intramolecular correlations (that is, correlations between sites belonging to a given m+xule), and p is the density of sites in the system. The symbol indicates the Fourier transform of a function. The SSOZ equation reduces to the OZ equation when G(k) is equal to the identity matrix. Neither the total correlation function, h, nor the direct correlation function, c, is known a priori. Therefore, in order to solve for h and c, another relation between h and c, known as the closure relation, is needed to complement eq 2.4. For simple fluids, the exact form of the closure relation is given byi4 [1

+ haMa‘M(r)l =

-@UaMa’M’(r)

where the index M refers to the type of molecule in the system, the index a refers to a particular site in the molecule, r is the distance between sites aM and d M ’ , uaMdMr(r)is the interaction energy between sites aM and dM‘,zaMis the charge on site aM, and AaMatMt and CUM,,, are parameters of the model associated with the strength of the interactions contributing to the Len-

10-3Aa~a~I kcal AI2/mol

+ haMa’M’(r)- C a M u d r ) + EaMa’M’(r) (2.5)

where EaMa*M,(r)are the bridge functions, which are complex functionals of the total correlation functions, and B = l / ( k T ) where k is the Boltzmann constant and T i s the absolute temperature. In the case of simple fluids, approximations enter the theory when a form for the bridge functions, ea^&), is chosen. In general, the pair correlation function, g(r) = h(r) + 1, can be

Lue and Blankschtein

8584 The Journal of Physical Chemistry. Vol. 96, No. 21, 1992

representedt6 as an infinite sum of diagrams, and the bridge functions are an infinite sum of graphs which are a subset of these diagrams. The bridge diagrams are highly connected and therefore are of a shorter range than the other diagrams which contribute to the total correlation function^.'^ If the bridge diagrams are neglected, that is, if one sets E(r) = 0

(2.6)

in eq 2.5,then one obtains the hypemetted-chain (HNC) closure. The HNC closure has been found to work well for fluids with Long-ranged forces such as those described by the restricted primitive model.’I There have been several approximations to the bridge functions of the form E = E(t(r)),where t(r) = h(r) - c(r), which have proven successful for simple fluids. If the bridge functions are approximated by E(r) = In [ l + t(r)]- t(r) (2.7) then one obtains the Percus-Yevick (PY) closure, which has been found to work well for systems with only short-ranged forces such as the hard-sphere and Lennard-Jones systems.” Note that the PY closure sums over fewer diagrams than the HNC ~ l o s u r e . ’ ~ The Martynov-Sarkisov (MaS) closure was proposed3as an improvement to the PY clamre. It sums all the diagrams contained in the HNC closure and approximates the bridge functions by E(r) = [I

when all the sites in a molecule are collapsed to a single point.2E In spite of these limitations, theories which employ the SSOZ equation with a closure of the form given in eq 2.5 have been found to perform well for a wide variety of moledar system~.~+’JS~~JC’J~ These theories have been found to yield pair correlation functions in good agreement with the results of computer simulations. In fact, the SSOZ equation is the most extensively used approach to model interaction-site fluids. Accordingly, in this paper, we will utilize this theoretical approach to study hydrocarbon-water systems. C. C d c u l a t h of Bulk Tbermudynamic Properties from tbe Pair Correhbloa Fuoctio~m.In principle, all the bulk thermodynamic properties of a system can be calculated once the various site-site pair correlation functions, g(r) = h(r) 1, are known. In Appendix B, we present mathematical expressions relating various useful bulk thermodynamic properties to the sittsite pair correlation functions, which will be utilized in section 3 and Appendix A to predict the thermodynamic behavior of the hydrocarbon-water and pure water systems, respectively. Here, we will focus on the development of a new analytical expression for the residual chemical potential, appropriate for interaction-site fluids, in terms of pair and direct correlation functions at full coupling for the various closure relations examined in this paper. The residual chemical potential of a molecule of type M, p g , is given by65

+

+ 2t(r)]t/2- t ( r ) - 1

(2.8) For the case of hard spheres, the MaS closure predictions of thermodynamicproperties, such as the pressure, as well as of the pair correlation functions, are in better agreement with computer simulations than the corresponding PY closure prediction^.^ The BallonePastore-Galli-Gazzillo (BPGG) closure was introduced4 as a generalization of the M a s closure and has an adjustable parameter, s. Its bridge function is given by E(r) = [I + ~t(r)]I/~- t(r) - 1 (2.9) It is interesting to point out that when s = 1, the BPGG closure reduces to the HNC closure, and when s = 2 it reduces to the MaS closure. In previous appIications,4s was fixed by enforcing the equality of the pressures obtained from the virial and compressibility equations. The Rogers-YoungI8 (RY) closure has an adjustable function, fir). Its bridge function is given by

The only restrictions on the function f i r ) are that it should satisfy the boundary conditions: fi0)= 0 andf(m) = 1. Thus, as can be seen by comparing eq 2.10 with eqs 2.7 and 2.6 in the limits fi0)= 0 andfi-) = 1, respectively, the RY closure mimics the PY closure at short distances, while it mimics the HNC closure at long distances. Typically, the functionflr) is chosen to be of the formfir) = 1 - exp(-ar). The parameter a is chosen such that the pressures calculated using the virial and compressibility equations are equal. The RY closure has been found to work much beter than the PY closure for systems with purely repulsive interactions, such as l/r“potential fluids.’* All the closure relations presented above were developed for simplefluids. Hirata and Rossky suggestedi9using the closure relations for simple fluids, in the context of the SSOZ equation, to describe the interaction-sitefluids. Note, however, that in the case of polyatomic molecules, using closures of the form given in eq 2.5 is in itself an approximation. Indeed, this will generate20 an infinite number of diagrams which are unallowed in the exact theoryZ1sz2 of interaction-site fluids at every level of the density.23 Approximations of this type lead to a theory with well-known and well-documented deficiencies. For example, (i) the predicted dielectric constant is always equal to that of an ideal gas,24,25 (ii) the total correlation functions depend on the presence of auxiliary ~ites,2~-~’ and (iii) the theory does not reduce to the proper limit

r

where couples the interaction of a single molecule of type M with the rest of the molecules in the system, a runs over the interaction sites of a molecule of type M, M’ runs over all molecular types, a’ runs over the interaction sites of molecules of is the interaction energy of a single molecule type M’, and ( of type M with the rest of the molecules in the system. Note that at infinite dilution, the solute-solvent interaction energy,, :e is given by

€,“v = ua‘M’pu.Mt$4r9 dr UaMa’M’(r)&Ma’M’(r) (2.12) where a runs over all the interaction sites of the solute molecule (which is of type M), M’ runs over all molecular species in the system, and d runs over the interaction sites of molecules of type M’. The only restrictions on the integration path in eq 2.1 1 are that (1) U , M d M ’ ( r ; r = l ) = uaMdM*(r)and that (2)u,MdM’(r;r=o) e 0. Note that the integration with respect to ‘j is performed at constant T, pa^, and V. In principle, one can use the expression in eq 2.1 1 to calculate the residual chemical potential by (a) solving the SSOZ equation for different values of the coupling parameter, (b) calculating the solute-solvent interaction energy, ( via eq 2.12,and (c) performing the coupling constant integration numerically. Clearly, it would be much more efficient if one could find an expression for p g which is only a function of the correlation functions at = 1. Then, the SSOZ equations would only have to be solved once for = 1. To accomplish this goal, we have g e n e r a l i i the derivation of Morita and Hiroike,’2 which applies to simple fluids only, to the case of interaction-site fluids. If a closure relation of the form given in eq 2.5 is used, then generalizing the arguments of Morita and Hiroike, p C m be expressed as follows:33

r,

r

r

For the HNC closure, where E(r) = 0, eq 2.13 reduces to the

Liquid-State Theory of Hydrocarbon-Water Systems

The Journal of Physical Chemistry, Vol. 96, No. 21, 1992 8585

expression obtained34by Singer and Chandler for interaction-site fluids. For other closures, however, the last term in eq 2.13, which requires the integration of the bridge functions over the coupling constant poses a computational problem. Note that for an approximate closure, the residual chemical potential is, in general, dependent on the path of integration chosen to evaluate the last term in eq 2.13. For one-component simple fluids, Kjellander and Sarman were able to find35an analytical expression for this integral for a class of closures of the form E = E(t(r)),where t(r) = h(r)- c(r),using a path along which h(r,r) = rh(r). Kjellander and Sarman’s evaluation of the last term in eq 2.13 can be extended directly to interaction-site fluids (see Appendix C) to yield the following result:”

r,

0.‘51

0.0

,

0

2

i

I

,

4

6

L,

6

r (4

Substituting eq 2.14 into eq 2.13, we obtain an analytical expression for the residual chemical potential in terms of the quantities evaluated at = 1 only. We believe that the lack of an analytical expression for the residual chemical potential (of the type given in eqs 2.13 and 2.14), for closures other than the HNC closure, may have discouraged, so far, the use of other closures to study the thermodynamic behavior of molecular fluids. The residual chemical potential at infinite dilution is directly related to the Henry’s law constant, KM, as follows:

1

r

KM = PkT exp[fl~Fl

(2.15)

where p is the density of pure solvent at the system temperature and pressure. The Henry’s law constant, which is an experimentally measurable quantity, is related to the solubility of gases in liquids by the following equation: KM =

PM~M

-

(2.16)

80.6

0’41

I 10

0.2

0.0

_2

1

,

I

4

8

6

10

r (4

Figure 1. Methane-water pair correlation functions (a, CH4-0, b, CH4-H) as a function of r for methane at infinite dilution in water at 25 OC and 0.997 g/c”: (i) SSOZ equation with HNC closure (-), (ii) SSOZ equation with MaS closure (-- -), and (iii) SSOZ equation with BPGG closure, s = 1.79 (.-).

’YMXM

aM

where p Mis the partial pressure of the solute (species M), is the fugacity coefficient of the solute in the gas phase, y M is the solute activity coefficient in the liquid phase, and XM is the solute mole fraction in the liquid phase. The solute activity coefficient, YM, is chosen such that it will tend to unity as XM approaches zero. Because the solubility of hydrocarbons in water is extremely low, yMshould be close to unity. Furthermore, at the conditions of interest in this work (p = 1 atm, T = 5-80 fC), the gas phase can be treated as being ideal.69 Accordingly, 4~ can be set equal to unity. Incorporating these approximations into eq 2.16, we find that the Henry’s law constant, KM, provides a simple relation between the solubility of a gas in a liquid, xM, and its partial pressure, p ~ Specifically, . KM = P M / X M . In section 3 and A p pendix A, we present the results of predictions of p$ for hydrocarbons at infinite dilution in water and pure water, respectively, using eq 2.1 1 or eqs 2.13 and 2.14. In principle, although not reported in this paper, the predicted values of p@ can be utilized to predict the Henry’s law constant, K M ,utilizing eq 2.15 and the solubility, xM,utilizing the relation xM = pM/KM. D. Computational Details. The SSOZ equations with a particular closure relation form a set of N(N 1)/2 coupled, nonlinear integral equations, where N is the number of distinct interaction sites. Note that both water and propane have two distinct interaction sites, and both methane and ethane have one interaction site. These equations were solved using an extension of the method developed by Labik et al.36 The grid spacing was 0.018 A with 4096 grid points. Halving the number of grid points changed the computed thermodynamic properties by less than 1%. The long-ranged interactions were treated using the method of Ng.37 The integration over the coupling constant to obtain p? in eq 2.1 1 was performed numerically, using an eight-point GaussLegendre quadrature. Doubling the number of quadrature points altered the value of the integral by less than 1%. The derivative of pC with respect to temperature in eq B7 was computed using a central difference approximation with a AT

+

value of 1 OC. Halving the AT value produced no significant change in the values of the computed derivatives. 3. Results

A. Pure Wnter. We solved the SSOZ equations with the HNC closure for the SPC and other potential models of water at ambient conditions ( T = 25 OC and p = 0.997 g/cm3); see Appendix A for complete details. We utilized the HNC closure because it has been found to work well for simple fluids with long-ranged interactions. B. M e h n e in Wnter. The SSOZ equations were solved for methane at infinite dilution in water with the HNC, PY, RY, MaS, and BPGG closures for the methanewater pair correlation functions and the HNC closure for the water-water pair correlation functions. We used the OPLS for methane and the SPC potential model for water. The CH4-O pair correlation functions are shown in Figure la, and the CH4-H pair correlation functions are shown in Figure 1 b. The solid lines are the predictions of the SSOZ equation with the HNC closure, the dashed lines are the predictions with the M a s closure, and the dotted lines are the predictions with the BPGG closure for s = 1.79. This s value was chosen such that the predicted p$ fits the experimental value at T = 25 OC and p = 0.997 g/cm3. Note that we have found that the SSOZ equation with the RY and PY ciosures predicts a negative region in the CH4-H pair correlation function. This is unphysical and reflects the well-known limitations of these closures when describing systems exhibiting strong attractive interactions. Accordingly, we do not report the results for these two closures. We were unable to find any methane-water pair correlation functions from simulations that utilized potentials similar to ours. In addition, to the best of our knowledge, no scattering experiments have been performed on methane-water systems. Consequently, no comparison with our theoretical predictions could be made. Predictions for the enthalpy of solution, ARM (see eq B7), the residual chemical potential, p g (see eq 2.1 l), the solute partial

Lue and Blankschtein

8586 The Journal of Physical Chemistry, Vol. 96, No. 21, 1992

TABLE II: OPLS Metbaw at Infdte Dilution in SPC Water at T = 25 OC red p = 0.997 p / c d A&,

kcal/mol -3.30," -3.15' -5.02d 1.56 -2.27 -1.98

exptl comp simul HNC MaS BPGG (s = 1.79)

kcal/mol 2.01' 2.3e 8.72 (8.71) 1.21 (4.64) 2.02 (5.15)

M

VhIl

CUV,

A)

kcal/mol

Nbyd

57.3; 61MC 53.6 65.2 63.1

-2.89d -1.75 -2.60 -2.50

19,' 20.3d 19.0 18.7 18.7

"See ref 38. bSee ref 39. CSeeref 40. dSee ref 42. eSee ref 43. 'See ref 41. 'See ref 44. 3.0 E

U

0.0

,

,

1

20

40

1

1

I

60

Bo

100

T ("C)

Figure 2. Residual chemical potential of methane at infinite dilution in water as a function of temperature: (i) predictions of the SSOZ equation with the MaS closure (---), (ii) predictions of the SSOZ equation with the BPGG closure, s = 1.79 (-), (iii) data of Morrison and Billett4 (O), (iv) data of Wen and Hung4s (w), and (v) data of Yaacobi and BenNaim4' (A).

molar volume, V,,(see eq B3), the solutedvent interaction energy, e,"v (see eq 2.12), and the hydration number, Nhydr of methane in water at T = 25 OC and p = 0.997 g/cm3 are listed in Table 11. Note that following Jorgensen et a1.,4*the hydration number, Nhydr of methane is taken to be the number of water oxygen sites within 5.35 A of the center of a methane molecule. The predicted values are compared with the results from experiments3841and computer simulation^.^^^ To demonstrate the path dependence of ~ 2we, calculated this quantity using two different integration paths: (1) h(r,f') = rh(r), and (2) u(r,r) = ru(r). The first entry in Table I1 for is the prediction using the first integration path, and the entry in parentheses is the result from the sacond integration path. In order to calculate p C using path 1, the SSOZ equations were solved at full coupling, and then eq 2.13 with eq 2.14 was applied. Because no such formula exists for path 2, the integral with respect to had to be performed numerically. This involves solving the SSOZ equation at each value of determining (e;),and then using numerical quadrature to compute the integral. Therefore, the calculations using path 1 were much less time consuming than those using path 2. As can be seen from Table 11, the HNC closure shows no dependence on the path of integration, while for the M a s and BPGG closures the first path yields a lower value for &' than the second path. It is noteworthy that the first path, which is the more computationally efficient one, is in better agreement with the experimental and computer simulation values. If the closure relation were exact, then &' would be independent of the integration path.35 However, because of the approximations in the theory, this dependency arises. This path dependence is analogous to the discrepancy which arises in the calculation of the pressure using the compressibility and the virial routes with an approximate closure. In Figure 2, we present the predicted temperature variation of the residual chemical potential of methane at infinite dilution in water over the range 5-80 OC. Again, we used the OPLS for methane and the SPC potential model for water. The HNC closure was used for the water-water pair correlation functions, and the variation of the density of water with temperature was taken from the data of Kell.45 The solid line is the prediction of the SSOZ equation with the BPGG closure with s = 1.79, and

pz

r

r,

, 0

1

6

4

8

10

r (A)

Figure 3. Ethantwater pair correlation functions (a, C H 4 ; b, CHI-H) as a function of r for ethane at infinite dilution in water at 25 O C and 0.997 g/cm3: (i) SSOZ equation with HNC closure (-), (ii) SSOZ equation with MaS closure (- -), and (iii) SSOZ equation with BPGG closure, s = 1.96 (-.). Note that the predictions of the MaS and B E G closures are almost identical.

-

the dashed line is the prediction of the SSOZ equation with the MaS closure. The value of the free parameter s of the BPGG closure was chosen to fit the experimental residual chemical potential at 25 OC and p = 0.997 g/cm3. The symbols are experimentally measured values.4648 As can be seen from Figure 2, the predictions of the BPGG closure are in very goad agreement with the experimental data. The MaS closure shows the same trend as the BPGG closure but is shifted downward. It is noteworthy that although not shown in the figure, the SSOZ-HNC equation grossly overpredicts (by more than 300%)the values of the residual chemical potential. The same calculationswere a h repeated for methane in water using the TIP3P model for water (see Table VII). The difference in the methantwater pair correlation functions using the SPC and TIP3P potential models of water is too small to be seen. In addition, we found that the choice of water potential model does not have a significant effect on the predictions of the bulk thermodynamic properties listed in Table 11. C. Ethane in Water. We have also solved the SSOZ equation with the HNC, PY, RY,Mas, and BPGG closures for ethane at infinite dilution in water using the OPLS to describe the ethane interactions and the SPC potential model to describe the water interactions. The CH,-O pair correlation functions are shown in Figure 3a, and the CH3-H pair correlation functions are shown

The Journal of Physical Chemisrry, Vol. 96, No. 21, 1992 8587

Liquid-State Theory of Hydrocarbon-Water Systems 3.0

TABLE IIk OPLS Ethane at M i t e Dilution in SPC Water at T = 25 OC .ad 0

0.997

Q / C d

exptl -4.72,' -4.61d 1.84' comp simul -2oc HNC 1.33 12.3 (12.3) Mas -4.19 1.62 (6.96) 1.83 (7.08) BPGG (S = 1.96) -4.02

81.86 -4.8oC

69.8 82.0 81.6

-3.49 -4.49 -4.41

"See ref 38. bSee ref 40. 'See ref 42. dSee ref 41.

L !0.0[ , 0

,

,

Bo

40

20

,

,

,

Bo

the results available from experimental m e a s ~ r e m e n t sand ~~*~~~~ computer The predicted temperature variation of the residual chemical potential of ethane at infinite dilution in water over the range 5-80 OC is presented in Figure 4. As with methane in water, we find that the prediction of the BPGG closure (with s = 1.96) is in very good agreement with the experimental data for temperatures below about 50 OC. However, at the higher temperatures (T>50 "C), the agreement begins to worsen. The dashed line is the prediction of the SSOZ equation with the MaS closure, and as in the case of methane, it is displaced downward from that of the BPGG closure. As with methane, the SSOZ-HNC equation grossly overpredicts (by more than 60096) the values of the residual chemical potential. D. Propane In Water. The SSOZ equations with the HNC, PY, RY, Mas, and BPGG closures were also solved to predict the propane-water pair correlation functions. The propane interactions were modeled using the OPLS, and the water interactions were modeled using the SPC potential. The CH2-0, CH3-0, CH,-H, and CH3-H pair melation functions are shown in parts a-d of Figure 5 . The solid lines are the predictions of the HNC closure, the dashed lines are the predictions of the MaS closure, and the dotted lines are the predictions of the BPGG closure with s = 2.05. As before, the parameter s was chosen such that the predicted residual chemical potential is equal to the

I

100

T ("C)

Epprrr 4. Reaidual chemical potential of ethane at infinite dilution in water as a functionof temperature: (i) predictions of the SSOZ equation with the MaS closure (- -), (ii) predictions of the SSOZ equation with the BPGG c h r c , s = 1.96 (-), (iii) data of Morrison and Bill& (O), (iv) data of Wen and Hung" (m), and (iv) data of Yaacobi and BenNaim" (A).

-

in Figure 3b. The solid lines are the predictions of the SSOZHNC equations, the dashed lines are the predictions of the SSOZ-MaS equations, and the dotted lines are the predictions of the SSOZ-BPGG equations with s = 1.96. The parameter s was chosen such that the predicted residual chemical potential is equal to the experimentally measured value at ambient conditions. Note that the predictions of the MaS closure (dashed lines) are almost identical to those of the BFGG closure (dotted lines). As for methane, the RY and PY closures yield negative regions in the CH3-H pair correlation functions (for the same reason explained earlier), and therefore the results for these closures are not presented. We could not find pair correlation functions from either simulations or scattering experiments to compare with our predictions. Various predicted thermodynamic properties for ethane in water are listed in Table 111, along with

1.4 1.2

-

L

1.5

1

1.0 -

h

jO.8

dl 0.6 . 0.4 -

0.4u ,

,

0.2 0.0

-__2

2

4

6

8

10

r (4 Figure 5. Propanewater pair correlation functions (a, CH& b, CH,-O; c, CH2-H; d, CH3-H) as a function of r for propane at infinite dilution in water at 25 OC and 0.997 g/cm3: (i) SSOZ equation with HNC closure (-), (ii) SSOZ equation with MaS closure (---), and (iii) SSOZ equation with BPGG closure, s = 2.05 (--). Note that the predictions of the Mas and BPGG closures are almost identical.

8588 The Journal of Physical Chemistry, Vol. 96, No. 21, I992 1.6

Lue and Blankschtein

1

I

1.4 1.2 I 1 .o

0.4

0.2

1

1

__0

2

4

6

r

8

0

10

20

80

40

(4

80

100

T ec,

Figure 6. Propane-water pair correlation functions (CH2-O) as a function of r for propane at infinite dilution in water at 25 OC and 0.997 g/cm3 using the potential model of Ichiye and Chandler:'" (i) SSOZ equation with the HNC closure (-), (ii) SSOZ quation with MaS closure (---),and (iii) molecular dynamics simulation of Ichiye and Chandler'"

Figure 7. Residual chemical potential of propane at infinite dilution in water as a function of temperature: (i) predictions of the SSOZ equation with the MaS closure (---), (ii) predictions of the SSOZ equation with the BPGG closure, s = 1.96(-), (iii) data of Morrison and Billett4 (a), and (iv) data of Wen and Hung4*(0).

(e-).

~

~

_

_

_

_

TABLE I V OPE3 Propane at Infinite Dilution in SPC Water at T = 25 O C and p = 0.997 g/cm3 M M ,

kcal/mol

WE,

kcal/mol

VM?

A'

exptl comp simul HNC MaS BPGG

(S

-5.38,"-5.56d 1.96" 1 17.46 1 3c 2.20 17.5 (17.4) 92.7 -5.41 2.35 (10.2) 107.3 = 2.05) -5.61 1.98 (10.0) 107.8

1.5

M e"",

1

kcal/mol S

-6.6W -4.46 -5.76 -5.79

I

f l . O ;-

S I I

0.5

I

1

'See ref 38. bSee ref 40. cSee ref 42. dSee ref 41. experimentally measured value a t ambient conditions. Note that the predictions of the M a s closure (dashed lines) are almost identical to those of the BPGG closure (dotted lines). The results of the RY and PY closures are not reported because, as explained earlier, they yield negative regions in the CH2-H and CH3-H pair correlation functions. Ichiye and Chandler solved1othe SSOZ equation with the HNC closure for propane at infinite dilution in water and also performed molecular dynamics simulations to compare their predictions. In their work, they u t i l i i I o potentials which are different from ours. Accordingly, to compare with their simulation work, we solved the SSOZ equation with the M a s closure for propane in water using the same potential models'O that they employed. We find that the predictions of the M a s and the HNC closures for all the propane-water pair correlation functions are very similar, with the M a s closure yielding a slightly better agreement with the molecular dynamics simulation data. As an illustration, in Figure 6 we show a comparison of the predicted CH2+ pair correlation functions with the molecular dynamics simulation data. The solid lines are the predictions of the HNC closure, the dashed lines are the predictions of the M a s closure, and the dotted lines are the molecular dynamics simulation results. The other propane-water g(r)'s show similar trends and are therefore not reported. Predictions for NE, and eg are listed in Table IV along with the results from experimental measurement^^^^^^^^ and computer simulations."2 Figure 7 shows the predicted temperature variation of the residual chemical potential of propane at infinite dilution in water over the range 5-80 "C. The agreement between the predictions of the BPGG closure (with s = 2.05) and the experimental values is extremely good a t the lower temperatures ( T < 40 "C); however, it worsens at the higher temperatures ( T > 40 "C). Once again, the M a s prediction shows the same variation with temperatue as that of the B E G closure, although this time it is shifted upward. As with methane and ethane, the SSOZ-HNC equation grossly werpredicts (by more than 800%) the values of the residual chemical potential. The propane molecule has two nonequivalent interaction sites. We have analyzed the effect of the ESA by solving the SSOZ equation with the M a s closure for three different scenarios. The

mM, vM,

i - -.-..... 00 0

1 2

I 6

4

r

8

10

(4

Figure 8. Application of the equivalent-site approximation (ESA) to propane in water. (a) Propantwater oxygen pair correlation functions: (i) CH2-O,full potential (thin dashed line); (ii) CH3-0, full potential (thick dashed line); (iii) CH2-0, modified potential (thin dotted line); (iv) CH3-0, modified potential (thick dotted line); and (v) CH,-O, equivalent-site approximation (solid line). (b) Propane-water hydrogen pair correlation functions:(i) CH2-H, full potential (thin dashed line); (ii) CH,-H, full potential (thick dashed line); (iii) CH2-H, modified potential (thin dotted line); (iv) CHI-H, modified potential (thick dotted line); and (v) CH,-H, equivalent-siteapproximation (solid line).

first scenario treats the propane molecule as having two distinct sites using the OPLS interaction parameters. The second scenario treats a propane molecule as having two distinct sites, but with the CH, groups treated as CHI groups. The third scenario employs the ESA, which treats the molecule as having sites which interact as CHI groups and are all equivalent. As before, the OPLS were used to model propane, and the SPC potential was employed to model water. The sitesite pair correlation functions for each of these systems are shown in parts a and b of Figure 8. In each figure,the dashed

_

~

Liquid-State Theory of Hydrocarbon-Water Systems

The Journal of Physical Chemistry, Vol. 96, No. 21, 1992 8589

TABLE V: Ropne h SPC Water (MIS Closure) at T = 25 O C and = 0.997 g/cm3: A Test of the Equivalent-Site Approximation (=A)

p

scenario 1 scenario 2 scenario 3 (ESA)

rB, kcal/mol

Aj

2.35 3.32 3.42

107.3 108.1 104.9

3.0

1

I

kcaqmol -5.76 -4.64 -4.69

lines represent the sitwite pair correlation functions for the first s~enario,the dotted lines are for the second scenario, and the solid lines are the results of the ESA (the third scenario). The thin (dashed and dotted) lines are for the CHz-water correlations, and the thick (dashed and dotted) lines are for the CH3-water correlations. The predictions for &’, and z,”y for each of these scenarios are compared in Table V.

vM,

4. Discussion The simulation data listed in Tables 11-IV are a compilation of results from several researchers using slightly different interaction potential models and slightly different solution conditions. There are also some difficulties in obtaining certain thermodynamic quantities from computer simulations. For example, the evaluation of the enthalpy of solution requires4zthe performance of extremely long simulations in order to get the proper statistics. This is the main cause of the poor agreement between the simulation and experimental values for AllM. All these issues should be kept in mind when comparing the simulation results with the predictions of the SSOZ equations. For each of the hydrocarbons in water that we have studied, the HNC closure yields a gross overprediction of the solute residual chemical potential, b c , and the wrong sign for the enthalpy of solution, ARM.Furthermore, this closure also underpredicts the magnitudes of the solute partial molar volume, and the solute-solvent interaction energy, . :e The magnitude of this underprediction increases as the hydrocarbon chain length increases. In this respect, it should be emphasized that the HNC closure is the only one which has been used so far to predict the bulk thermodynamic properties of interaction-site fluids. This is probably due to the fact that until the present work, this has been the only closure for which an analytical expression for the residual chemical potential was available. The PY and RY closures both yield negative regions in the hydrocarbon-water hydrogen pair correlation functions for all the hydrocarbon-water systems that we have studied. As explained earlier, this reflects the well-known deficiencies of these closures when describing systems exhibiting strong attractive interactions. In addition, this may also reflect the overestimation3 of the contribution from the bridge functions to the pair correlation functions. On the other hand, the MaS and BPGG closures do a reasonable job in predicting the structural and bulk thermodynamic properties of hydrocarbon-water systems. We find that as the length of the hydrocarbon chain increases, the value of s required to fit the experimental p C increases. For methane, s = 1.79, for ethane, s = 1.96, and for propane, s = 2.05. Note that in each case, the values of s are closer to s = 2, the value for the MaS closure, than to s = 1, the HNC closure value. Recall from section 2 that the approximations made with the SSOZ-HNC equations lead to the inclusion of diagrams which are not present in the rigorous theory. In addition, the HNC closure neglects a class of diagrams, known as the bridge diagrams, which are present in the rigorous theory. Because the hydrocarbon-water interactions are short-ranged, the bridge diagrams, which are also short-ranged, are important. The MaS closure offers a better approximation for the bridge diagrams than the HNC closure and, consequently, yields better results. For methane and ethane in water, the prediction of the BPGG closure with the fitted parameter, s, is in very good agreement with the experimental data for the variation of the residual chemical potential with temperature over the relatively broad range 5-80 OC. Note that neither the interaction potentials nor the form

t

OS

I

0.0

1

I1

1

1

2

6

4

8

10

r (A)

i 9. Oxygen-oxygen pair correlation functions,goo(r), as a function of the sitesite separation,r, for SPC water at 25 OC and 0.997 gjcm’: (i) predicted using the SSOZ equation with the HNC closure (-), (ii) predicted from Monte Carlo simulations of Jorgensen et al.51(- - -), and (iii) deduced from neutron-scatteringmeasurements of Sopcr and PhilF

lip@ (-). 1.4

:I

b 0.8 -

d 0.6 1 i 0.4

vM,

0.0

’

2

4

1

I

I

6

E

10

r (A)

Figure 10. Hydrogen-hydrogen pair correlation functions, gHH(r), as a function of r for SPC water at 25 OC and 0.997 g/cm3. The notation is the same as in Figure 9. 2.c

1.5

v

51.0

0

0.5

0.0

0

2

6

4

8

I

r (4 Figure 11. Oxygen-hydrogen pair correlation functions, go&), as a function of r for SPC water at 25 OC and 0.997 g/cm3. The notation is the same as in Figure 9.

of the closure was adjusted with temperature. This agreement is somewhat surprising, considering the fact that the water-water pair correlation functions were only in qualitative agreement with the neutron scattering data (see Appendix A and Figures 9-11). This might suggest that only the short-ranged correlations in water, that is, the local packing of the water molecules around the hydrocarbon, are important in the hydrocarbon-water systems. For propane in water, the results for p r vs Tare still fairly good but not as impressive as for methane and ethane. The agreement breaks down at the higher temperatures ( T > 40 “C). The fact that the agreement worsens at the higher temperatures may reflect the inadequaciesof the SPC water model at these conditions. In

Lue and Blankschtein

8590 The Journal of Physical Chemistry, Vol. 96, No. 21, 1992 TABLE VI: Characteristic Properties of ma, A h a , deg 104.523" exptl 0.957 18" 109.47 SPC' 1 104.5 MCYd 0.9572 109.47 SPC/Ee 1 104.52 TIP3Pf 0.9572

the Water Molecule d, D ami, A' 1.856 1.4444 2.274 0 0 2.19 (2.438)' 0 2.351 0 2.347

"See ref 52. bSee ref 53. 'See ref 13. dSee refs 49 and 7. 'See ref 50. 'See ref 51. 8Value for the fitted 12-6-1 potential 3-site model. See text.

this respect, it is noteworthy that the deviation between the calculated vapor-liquid coexistence curve of water, obtained by computer simulation using the SPC potential model for water, and the corresponding experimental coexistence curve was founda to increase as the temperature increases. Figure 8 indicates that the propane-water pair correlation functions do not change significantly if one replaces the two CH, groups in the molecule by CHI groups. In addition, Figure 8 reveals that the propane-water pair correlation functions predicted using the ESA lie between the CH,-water and CH,-water pair correlation functions. This is reasonable, since one expects that these pair correlation functions should serve as some sort of average of the other two g(r)'s. For propane, the difference in the predictions of scenarios 1 and 2 (in Table V) is due solely to the difference in the interaction potentials. We see that this has a fairly significant impact. However, this difference will diminish as the hydrocarbon chain length increases. The discrepancy between the predictions of scenarios 2 and 3 (in Table V) is due to the inaccuracies of the equivalent-siteapproximation. According to Table V, the ESA does not have a significant effect on the prediction of the bulk thermodynamic properties. Specifically, the ESA predictions are within 5% of those obtained using the rigorous calculations. Moreover, it should be kept in mind that the error inherent in the ESA is most severe for propane. For ethane, the ESA is exact, and for the longer hydrocarbon chains, the ESA should become better. 5. Concluding Remarks Previous applications of the SSOZ equation to predict bulk thermodynamic properties of molecular fluids have employed

almost exclusively the HNC closure. In this paper, we have studied hydrocarbons (methane, ethane, and propane) at infinite dilution in water and have found that using the HNC closure for the hydrocarbon-water pair correlation functions yields poor results for the predictions of bulk thermodynamic properties. In addition, the PY and RY closures yield negative regions in the hydrocarbon-water hydrogen pair correlations functions, which is unphysical and reflects, as explained earlier, the inherent limitations of these closures when describing systems exhibiting strong attractive interactions. On the other hand, using the MaS or the BPGG closures, which have not been used previously with the SSOZ equation, for the hydrocarbon-water pair correlation functions yields fairly good predictions of the structural and bulk thermodynamic properties. In fact, the temperature dependence of the residual chemical potential over the relatively broad range 5-80 OC is captured very well. We have also extended the work of Kjellander and S a r ~ n a n , ~ ~ aimed at calculating residual chemical potentials, to interaction-site fluids. For approximate closures, the residual chemical potential is, in general, dependent on the integration path. Along the path, h ( r ; r ) = r h ( r ) , we were able to find an analytical expression for p r in terms of correlation functions at'j = 1, for closures of the form given in eq 2.5 and E = E(t(r)),where t(r) = h(r) - c(r). This new analytical expression facilitates the calculation of the residual chemical potential by eliminating the previous need to perform a numerical integration over the coupling constant. We hope that the availability of this new analytical expression for will encourage the use of other closures, in addition to the HNC closure, in future studies of the therm4ynamic behavior of molecular fluids. We have also tested the predictive accuracy of the equivalent-site approximation on propane and found that it performs fairly well. Note that of all the n-alkanes, propane ~OSCSthe most severe challenge to the ESA. The ESA was found to yield predictions of the bulk thermodynamic properties which are within 5% of those obtained using the rigorous calculations. We are currently investigating longer hydrocarbons in water in order to predict their structural and bulk thermodynamic properties, as well as to examine solution effects on the hydrocarbon chain conformations.

fir

TABLE MI: Interaction Parameters of the Potential Model8 of WateP ZH

SPCb MCY' SFC/Ed TIP3P

(-20/2) 0.41 0.433 04 0.423 8 0.417

A009

(kcal/mol) AI2 624 OOO 309 408 629 400 582 OOO

Aon, (kcal/mol) A'* 0 689.348 0 0

AHH (kcal/mol) A'* 0 610.455 0 0

COO,

(kcal/mol) A6 625.5 -262.566 625.5 595.0

OCOHand CHHare equal to zero. bSee ref 13. 'See ref 7. dSee ref 50. 'See ref 51. TABLE MII: Predictions (Bolded h t r y la h c b Cell) of the SSOZ-HNC Equrtions for W8ter at 25 "C and 0.997 g/cm3 d Comprrbm with Computer Simulrtion d Experimental Values expt MCY SPC SPC/E TIP3P -P, kcal/mol 9.9v 8.80 9.90 9.90 (11.1) 9.91 8.58: 8.51' 10.18,c 9.01,' 10.8' 9.8@ 9.82) 3.74 5.15 4.93 (6.18) 5.22 -Aw, kcal/mol 5.74" 4.31' 5.8e 5.5 (6.4)' 5.4' -p,kcal/mol 6.324b 2.03 2.50 2.52 (3.77) 2.70 1O6r, atm-' 45.2472c 30.3 55.6 57.1 58.9 53' 278 41.41h 188 Cy, -1 g-' K-' 0.9899 0.932 0.869 0.911 0.859 0.999 1.03,'0.909,' 1.14'

"See!ref 57. 6See ref 58. 'See ref 45. dSee ref 66. 'NVT Monte Carlo simulation with T = 25 "C and p = 1 g/cm3. See ref 59. 'NVT Monte Carlo simulation with T = 25 "C and p = 0.997 g/cm3. See ref 61. BNPT Monte Carlo simulation with p = 1 atm and T = 25 "C. See ref 51. *NPT molecular dynamics with p = 1 atm and T = 25 "C. See ref 62. 'See ref 63. 'See ref 64. For an explanation of the entries in parentheses for the SPC/E model, see Appendix A.

Liquid-State Theory of Hydrocarbon-Water Systems Acknowledgment. This research was supported in part by the National Science Foundation (NSF) Presidential Young Investigator (PYI) Award to Daniel Blankschtein and an NSF Grant (No. DMR-84-18718) administered by the Center for Materials Science and Engineering at MIT. Daniel Blankschtein is grateful for the support of the Texam-Mangelsdorf Career Development Professorship at MIT. He is also grateful to the following companies for providing PYI matching funds: BASF, British Petroleum America, Exxon, Kodak, and Unilever. Leo Lue is grateful for the award of an NSF Graduate Fellowship. Finally, we thank J. G. Harris for many illuminating discussions.

Appendix A: Predictions of Structural and Bulk Thermodynamic Properties of Pure Water 1. Potential Models of Water. We have used the SPC," MCY,49 SPC/E,SOand TIP3PS' potentials to model water. Some characte!istic Properties associated with these four potential models are listed in Table VI, along with the mmponding experimentally d e t e r ~ n i n e dvalues ~ ~ ~ ~for~ isolated water molecules. In Table VI, f O H is the oxygen-hydrogen bond length, OHOH is the hydrogenoxygen-hydrogen bond angle, d is the dipole moment, and cypl is the polarizability. The experimental values are for water molecules in vacuo. The parameters z, A, and C corresponding to these four potential models of water are listed in Table VII. Note that the original MCY model has49four interaction sites, that is, one on each of the three atoms in the water molecule and a fourth site located at the center of the molecule, directly below the oxygen atom. Furthermore, the sittsite interaction potentials are complicated functions of the sitesite separations. In view of these complexities, Pettitt and Rossky have fit7the MCY energy surface to a three-site model, thus obtaining a modified MCY model which is similar in mathematical structure to the other three potential models of water. This is the form of the MCY model which we have used. In each of the potential models of water, with the exception of the MCY model, there is no oxygen-hydrogen repulsion but only a Coulombic attraction (AoH = 0; see Table VII). This will lead to a "catastrophic overlap- of the oxygen and hydrogen sites of different water m~lecules.~ This phenomenon is not observed in simulation studies because there is a large energy barrier of about lo6 kcal/mol that must be overcome in order to bring together the oxygen and hydrogen sites of two different water molecules. The simulation is set up initially such that all the water molecules are located sufficiently far away from each other such that they are outside this "collapse barrier". Since the typical kinetic energy of a simulated water molecule at room temperature is about 1 kcal/mol, the region of phase space where the oxygen and hydrogen sites overlap has a low probability of being sampled within a limited-time simulation. Note, however, that if a simulation were started with the water molecules located sufficiently close to each other such that they were inside the energy barrier, then this collapse would be observed. In the case of the SSOZ equation, which in effect samples all phase space, this collapse of the water molecules would occur in the absence of an oxygen-hydrogen repulsion. To prevent this, we have added' a 1/rI2repulsive term having a magnitude, AoH, chosen so as to reproduce the experimentally observed value of P.Specifically, we find that for the SPC model, AOH = 900 AI2 kcal/mol, for the SPC/E model, AoH = 800 AI2 kcal/mol, and for the TIP3P model, AoH = 1030 AI2kcal/mol. These values are of the same order of magnitude as the value chosen7by Pettitt and Rossky (AOH= 226 AI2 kcal/mol). 2. S t " d Rojwrtics. We have solved the SSOZ equations with the HNC closure and soft mean-spherical approximation (SMSA) closure (see eq A2) for the SPC, MCY, SPC/E, and TIP3P potential models of water (see Tables VI and VII) at ambient conditions ( T = 25 "C and p = 0.997 g/cm3). We have employed these closures because they have been found to work well for simple fluids with long-ranged interactions. The soft mean-spherical approximations4(SMSA) closure is not of the form given in eq 2.5. Instead, it requires that the potential function, u(r), be divided into a repulsive contribution,

The Journal of Physical Chemistry, Vol. 96, No. 21, 1992 8591

uo(r),and an attractive contribution, ul(r);that is, u(r) = U O ( 4 + U l ( 4 (All This division is typically done using the Weeks-Chandler-AndersenSS(WCA) criterion. The SMSA closure is then given byS4

+

1 + h(r) = e-8uo(r)[l h(r) - c(r) - &(r)]

(A2)

The SMSA closure has been found to work extremely well for Lennard-Jones fluids and provides a significant improvement over the HNC and PY closures for liquid metals.S4 Plots of the sittsite pair correlation functions, goo(r), gHH(f), and go&), for the SPC potential model of water are given in Figures 9-1 1. The solid lines are the predictions of the SSOZ equation with the HNC closure. The dashed lines are the results from Monte Carlo simulations performed by Jorgensen et aL5' The dotted lines are the g(r)'s deduced from neutron-scattering measurements performed by Soper and Phillips.56 We have found that the SMSA closure predicts a negative region in both the hydrogen-hydrogen and oxygen-hydrogen pair correlation functions. This is unphysical, and therefore, we do not report the predictions based on this closure. We have also obtained the pair correlation functions for the MCY, SPC/E, and TIP3P models for water. The g(r)'s for these water models were found to be very similar to those for the SPC model and, therefore, are not reported. The agreement between the goo(r) predicted by the SSOZHNC equations (full line) and that predicted by the simulation data (dashed line) is rather poor; see Figure 9. The position of the first peak is approximately correct, but the positions of the other peaks are shifted. There is an indication of a peak in the predictions of the SSOZ-HNC equations at the location of the second peak in the experimental and simulated goo(r)'s. Overall, the predicted structure is reminiscent to that of a Lennard-Jones fluid, a hexagonally closed-packed structure. Recall, however, that water is actually tetrahedral in structure. The agreement of the predicted gHH(f) and goH(f) with the simulation and experimental results is fair; see Figures 10 and 11. The peaks which are shown for the H H and OH g(r)'s are due to correlations between sites located on adjacent water molecules; that is, they reflect interactions between a water molecule and its first coordination shell. It appears, therefore, that correlations between water molecules at short range are captured fairly well by the SSOZ-HNC equations. The incorrect prediction of the "long-ranged structure" (see Figure 9) is possibly due to the known defects of the SSOZ equation, in particular, the inability to predict nontrivial angular correlation factors. These angular correlations are important in water because of its highly directional hydrogen-bonded structure. One way to improve the liquid state theory predictions may be to use the Chandler, Silbey, and Ladanyi (CSL)equati~ns,~ which do yield nontrivial angular correlations. 3. Bulk Thermodynamic Properties. Table VI11 lists the predicted values of the residual internal energy, P (see eq Bl), the residual Helmholtz free energy, A" (see eq BS), the residual chemical potential, p z (see eq 2.1 l), the isothermal compressibility, K (see eq B2), and the heat capacity at constant volume, C y (see eq B4), for water at 25 "C and 0.997 g/cm3 using the SSOZ-HNC equations for the various potential models of water. Also listed are the values of these properties obtained from experimental measurement^^^^^^^^^^^ and computer simulation^.^^ In Table VIII, the first column contains the experimentally measured values of each property. Each of the four remaining columns corresponds to predictions of these properties for various potential models of water. For each property, the bolded entry in each cell corresponds to our predictions based on the SSOZHNC equations, and the remaining entries correspond to reported calculations from various computer simulations. In the four potential models for water described above, the molecules cannot alter their dipole moments and consequently uN = 0 (see Table VI). From Table VI, we also see that each of the dipole moments of the model water molecules is larger than the actual (experimental) dipole moment of an isolated water molecule.

Lue and Blankschtein

8592 The Journal of Physical Chemistry, Vol. 96, No. 21, 1992

For the SPC/E model, Berendsen et al. attributeSothis to the polarization of the molecule by its environment. Using simple arguments, they find that the energy associated with polarizing a water molecule, EPI,is given byS0

E,I =

reading section 3, in this Appendix we present expressions for each of the thermodynamic properties which we have calculated. The residual internal energy, P, is given by65

( d - do)*

2ffPl

where d is the dipole moment of a model water molecule, d’ is the dipole moment of an isolated water molecule, and ap0lis its polarizability, both taken to be equal to their experimental values in vacuo. This offers a method to correct for the polarizability of the SPC/E water molecules. In Table VIII, the bolded entries in parentheses for the SPC/E model correspond to values predicted without the polarization correction, As explained earlier, in the case of the SPC,SPC/E, and TIP3P water potential models, a repulsive term was added’ to prevent the overlap of the oxygen and hydrogen sites. The magnitude of this repulsive term was estimated by demanding that the predicted value of P be equal to the experimentally measured one. This is clearly represented in the bolded entry values of W e (-9.90 kcal/mol) for these three water models. Note, however, that the prediction of P for the MCY water model (-8.80 kcal/mol) is an independent one. The major difficulty in simulating water is in handling the long-ranged interactions. The variation in the reported simulation values presented in Table VI11 is due mainly to the different ways in which these long-ranged interactions are treated. For example, the values -9.01 and -10.8 kcal/mol for P of SPC water are a result of employing two different smoothing functions to truncate the interaction ~ t e n t i a l . 6Using ~ Ewald summation yields,@P = -9.82 kcal/mol. This variation in the simulation values should be kept in mind when comparing the results of the simulations with the predictions of the SSOZ equation. Overall, an examination of Table VI11 indicates a fair agreement between the predictions of the SSOZ-HNC equations and the experimental results, with the exception of The fact that the prediction for Are is good and that = A / N p for a one-component system (where p is the chemical potential, A is the Helmholtz free energy, N is the number of moles, and p is the pressure) suggests that the poor prediction of the residual chemical potential may be related to the fact that pressure is very sensitive to the approximations made in liquid-state theory. The residual Helmholtz free energy was computed using eq B6, as well as by numerical integration as spedied by eq BS. By using an eight-point Gauss-Laguerre quadrature, the two methods yielded predictions within 1% of each other. This is a reflection of the accuracy of our numerical methods. Mezei et al. computedsg A“ for MCY water by performing a Monte Carlo simulation at each value of the coupling parameter, {, to obtain ( Va)p and subsequently utilizing an eight-point Gauss-Legendre quadrature to perform the integration with respect to the coupling constant (see eq BS). When we used an eight-point Gauss-Legendre quadrature, we obtained Am = -4.46 kcal/mol, close to their value of -4.31 kcal/mol. However, we found that if we increased the number of quadrature points, the computed value of Are slowly increased to -3.74 kcal/mol, the value reported in Table VIII. Therefore, we feel that the results of Mezei et al. would be closer to ours if they had used more quadrature points. We seemed to have captured the short-ranged correlations in water properly. However, the distribution of water molecules in the second coordination shell and beyond, as predicted by the SSOZ-HNC equations, was found to deviate significantly from that found experimentally (see Figure 9). This may not pose a problem because the water-hydrocarbon interactions are shortranged, and therefore, the hydrocarbon is not expected to interact significantly with water molecules beyond the first coordination shell.

+

Appendix B Formulas for tbe T h e r m o d y ~ mProperties i~ As stated in section 2, once all the site-site pair correlation functions are known, all the thermodynamic properties of the system can be computed. For completeness and to facilitate

where Vis the system volume, p is the density of sites in the system, M and MI run over all molecular species, and a! and a!/ run over all the sites in each molecule. Here, the residual of a property is defined as the difference in the value of a particular system property from that corresponding to an ideal gas at the same temperature, composition, and total volume. Note that this definition of the residual is slightly different from that typically used in classical thermodynamics, where the residual is defined with respect to an ideal gas at the same temperature, mole numbers, and pressure.66 For a one-component system (for example, pure water), the isothermal compressibility,i,is given by6’ K

= ‘[I P

+ pi(())]

(B2)

where p is the system density. For a two-component system, the solute partial molar volume, Vu,is given by6’

where pu is the solute density and pv is the solvent density. The index u refers to the solute, and the index v refers to the solvent. The heat capacity at constant volume, Cy,is given by66

where C$,$is the ideal-gas heat capacity for molFules of type M, which depends on tempemture only. Values of C$,$ for various molecules can be found in the The residual Helmholtz free energy Are, can be calculated as follows:6*

where the parameter {couples the interactions of all the molecules with one another and ( P)f is the residual internal energy of the system at a coupling of strength {. When { = 1, all the molecules in the system interact fully with one another; that is, U,M~M,(G{) = U,M!M’(r) for all a,a!’, M, and M’. On the other hand, when { = 0, there are no intermolecular interactions present in the system; that is, U,M~M~(I;{)= 0 for all a,d,M, and M’. In other words, in the limit { = 0, the system behaves like an ideal gas. Note that the integration with respect to {in eq B5 is performed at constant temperature, species densities, and total volume. For simple fluids, Morita and Hiroikesz were able to express Am in termsof integrals over the pair correlation functions. Singer and Chandler34generalized this result for the SSOZ equation under the HNC closure Generalizing both of these derivations for the SSOZ equation with a closure of the form given in eq 2.5, we find that eq B5 can be expressed as

1

s { l n (det[l - p+(k)Z(k)]) + Trp+(k)Z(k)J+

Liquid-State Theory of Hydrocarbon-Water Systems

The Journal of Physical Chemistry, Vol. 96, No. 21, 1992 8593

For the HNC closure, for which E(r) = 0, eq B6 yields an expression for ArQas a function of the correlation functions only at full coupling (t = 1). The enthalpy of solution of molecules of type M ARM,is given by69

where is the partial molar enthalpy of M in the ideal-gas state, PM is the partial molar enthalpy of M at infinite dilution, apis the coefficient of thermal expansion of the solution, and K~ is the isothermal compressibility of the solution.

Appendix C Derivation of Equation 2.14 If the bridge functions are approximated by the form E = E(r(r)), where t(r) = h(r) - c(r), then the last term on the right-hand side of eq 2.13 can be rewritten as

If we choose an integration path such that h ( r ; r ) = r h ( r ) for all solutesolvent pair correlation functions, we can write eq C1 as

(C2) Because the solute is at infinite dilution in the solvent, the solvent-solvent pair correlation functions are independent of the solutesolvent coupling parameters, Choosing h&r) = rhUv(r)implies that rUv(r;r) = rtuv(r)for all sites u associated with the solute and sites v associated with the solvent. This can be demonstrated by writing tuv(r)in terms of h(r) in the context of the SSOZ equation (eq 2.4); that is,

r.

?&r ) = iUVW - x ~ u ~ ’ ( k ) j i u v)[W) ( k r - Pd(k)l,.,-’ u’v‘

(C3) where u and u’ are sites associated with the solute and v and v’ are sites associated with the solvent. If the molecule is rigid, the only factor which has a dependence on is hut,,. Therefore, as stated before, we find that tuv(r;r)= rtUv(r)for all sites u associated with the solute and sites v associated with the solvent. We can now express as

r

r r = t Ut uh ~i uu~ ~~ 1M~ 4’ r )

Substituting eq C4 into eq C2, we obtain

This completes the proof of eq 2.14.

Reglatry No. Methane, 74-82-8; ethane, 74-84-0; propane, 74-98-6.

References and Notes (1) For comprehensive experimental and theoretical surveys of the hydrophobic effect, see.: (a) Ben-Naim, A. Water and Aqueour Solutions; Plenum: New York, 1974. (b) Tanford, C. The Hydrophobic Effect: Formation of Micelles and Biological Membranes; Wiley: New York, 1973. (c) Israelachvili, J. N. Infermolecular and Surface Forces; Academic: London, 1985; p 102. (2) Chandler, D.; Andersen, H. C. J . Chem. Phys. 1972, 57, 1930. (3) Martynov, G. A.; Sarkisov, G. N. Molec. Phys. 1984, 49, 1495. (4) Ballone, P.; Pastore, G.; Galli, G.; Gazzillo, D. Mol. Phys. 1986, 59, 275. ( 5 ) Curro, J. G.; Schweizer, K. S.J . Chem. Phys. 1987, 87, 1842. (6) Pratt, L. R.; Chandler, D. J . Chem. Phys. 1977, 67, 3683. (7) Pettitt, B. M.; Rossky, P. J. J . Chem. Phys. 1982, 77, 1451. (8) Yu, H.; Karplus, M. J . Chem. Phys. 1988,89, 2366. (9) Yu, H.; Row, B.;Karplus, M. J . Chem. Phys. 1990, 92, 5020. (10) Ichiye, T.;Chandler, D. J . Phys. Chem. 1988, 92, 5257. (1 1) Lee, L. L. Molecular Thermodynamics of Nonideal Fluids; Butterworths: Boston, 1988. (12) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J . Am. Chem. Soc. 1984, 106,6638. (13) Berendsen, H. J. C.; Postma, J. P. M.; von Gunsteren, W. F.; Her-

mans,J. Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981; p 331. (14) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic: London, 1986. (15) Chandler, D. J . Chem. Phys. 1973,59, 2742. (16) For an excellent review, see: Stell, G. The Equilibrium Theory of Classical Fluids; Frisch, H. L., Lebowitz, J. L., Eds.;Benjamin: New York, 1964; p 11-171. (17) Stell, G. Physica 1963, 29, 517. (18) Rogers, F. J.; Young, D. A. Phys. Reu. A 1984, 30, 999. (19) Hirata, F.;Rossky, P. J. Chem. Phys. Lett. 1991, 83, 329. (20) Chandler, D.; Silky, R.; Ladanyi, B. M. Molec. Phys. 1982,46,1335. (21) Ladanyi, B. M.; Chandler, D. J . Chem. Phys. 1975,62, 4308. (22) Chandler, D.; Pratt, L. R. J . Chem. Phys. 1976, 65, 2925. (23) Chandler, D. Molec. Phys. 1976, 31, 1213. (24) Sullivan, D. E.; Gray, C. G. Molec. Phys. 1981, 42, 443. (25) Chandler, D.; Joslin, C. G.; Deutch, J. M. Molec. Phys. 1982,47, 871. (26) Hsu, C. S.;Chandler, D.; Lowden, L. J. Chem. Phys. 1976,14,213. (27) Cummings, P. T.;Gray, C. G.; Sullivan, D. E. J . Phys. A 1981,14, 1483. (28) Chandler, D.; McCoy, J. D.; Singer, S . J. J . Chem. Phys. 1!W, 85, 5977. (29) Johnson, E.; Hazoume, R. P. J . Chem. Phys. 1979, 70, 1599. (30) Curro, J. G.; Schweizer, K. S.;Grest, G. S.;Kremer, K. J . Chem. Phys. 1989,91, 1357. (31) Hirata, F.; Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1982,77,509. (32) Morita, T.; Hiroike, K. Prog. Theo. Phys. 1960,23, 1003. (33) Note that a similar expression was obtained using the molecular

Ornstein-Zemike equation by: Chen, X.S.;Forstmann, F.; Kasch, M. J. Chem. Phys. 1991, 95, 2832. (34) Singer, S. J.; Chandler, D. Molec. Phys. 1985, 55, 621. (35) Kjellander, R.; Sarman. S . J. Chem. Phys. 1989, 90,2768. (36) Labik, S.;Malijevsky, A,; Vonka, P. Molec. PHys. 1985, 56, 709. (37) Ng, K. J. Chem. Phys. 1974,61, 2680. (38) Abraham, M. H. J. Chem. Soc., Faraday Trans. 1 1984,80, 153. (39) Reid, R. C.; Prausnitz, J. M.; Sherwood, T. K. The Properties of Gases and Liquids; McGraw-Hill: New York, 1977; p 358. (40)Moore, J. C.; Battino, R.; Rettich, T. R.; Handa, Y. P.; Wilhelm, E. J. Chem. Eng. Dafa 1982, 27, 22. (41) Dec, S. F.; Gill, S.J. J . Sol. Chem. 1984, 13, 27. (42) Jorgensen, W. L.; Gao, J.; Ravimohan, C. J. Phys. Chem. 1985,89, 3470. (43) Jorgensen, W. L.; Blake, J. F.; Buckner, J. K. Chem. Phys. 1989,129, 193. (44) Guillot, B.; Guissani, Y.; Bratos, S.J. Chem. Phys. 1943,95, 3643. (45) Kell, G. S.J. Chem. Eng. Data 1975, 20, 97. (46) Morrison, T. J.; Billett, F. J . Chem. Soc. 1952, 3819. (47) Yaacobi, M.; Ben-Naim, A. J . Phys. Chem. 1974, 78, 175. (48) Wen, W.; Hung, J. H. J. Phys. Chem. 1970, 74, 170. (49) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976,64, 1351. (50) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J . Phys. Chem. 1987, 91, 6269. (51) Jorgensen, W. L.;Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79, 926. (52) Eisenberg, D.; Kauzmann, W. The Strucfure and Properties of Wafer;Oxford: London, 1969. (53) Dyke, T.R.; Meunter, J. S . J . Chem. Phys. 1976,64, 1351. (54) Madden, W. G.; Rice, S.A. J . Chem. Phys. 1980, 72,4208. ( 5 5 ) Weeks, J. D.; Chandler, D.; Andenen, H. C. J . Chem. Phys. 1971, 54, 5237; Phys. Rev. A 1971, 4 , 1597. (56) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47. (57) Dashevsky, V. G.; Sarkisov, G. N. Molec. Phys. 1974, 27, 1271. (58) Ben-Naim, A,; Marcus, Y. J . Chem. Phys. 1984, 81, 2016.

8594

J . Phys. Chem. 1992, 96, 8594-8599

(59) Mezei, M.; Swaminathan, S.;Beveridge, D. L. J . Am. Chem. SOC. 1978, 100, 3255. (60) Hermans, J.; Pathiaseril, A.; Anderson, A. J . Am. Chem. Soc. 1988, 110, 5982. (61) Lie, C. G.; Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976, 64, 2314. (62) Motakabbir, K.A.; Bcrkowitz, M. J . Phys. Chem. 1990, 94, 8359. (63) Revost, M.; van Belle, D.; Lippens, G.; Wodak, S . Mol. Phys. 1990, 71, 587.

(64) de Pablo, J. J.; Prausnitz, J. M.; Strauch, H. J.; Cummings, P. T. J . Chem. Phys. 1990, 93,1355. (65) Hill, T. L. Statistical Mechanics: Principles and Selected Applicarions; McGraw-Hill: New York, 1956. (66) Smith, J. M.; Van Ness, H. C. Introduction to Chemical Engineering Thermodynamics; McGraw-Hill: New York, 1987. (67) Kirkwood, J. G.; Buff, F. P. J . Chem. Phys. 1951, 19, 774. (68) Hiroike, K. J . Phys. SOC.Jpn. 1960, IS, 771. (69) Wilhelm, E.; Battino, R.; Wilcock, R. J. Chem. Rev. 1977, 77, 219.

Effect of Restricted Geometries on the Structure and Thermodynamic Properties of Ice+ Y. Paul Handa,* Marek Zakrzewski? and Craig Fairbridges Institute for Environmental Chemistry, National Research Council of Canada, Ottawa, Ontario, Canada K I A OR6 (Received: May 18, 1992)

The thermal properties and structure of water in porous Vycor glass and various silica gels with pore radii in the range 23-70 A were investigated using calorimetry and X-ray diffraction. Samples containing pore water only and those containing pore water in equilibrium with bulk water were studied. The amounts of bound and freezable water in the samples were determined by measuring the heat of melting as a function of pore water content and by application of thermoporometry. The melting point of ice was depressed by 9-20 K, depending on the pore size. The heat of melting of pore ice was found to be 1648% smaller than that of bulk ice, depending on the pore size and whether the bulk phase was present or not. On quenching to liquid nitrogen temperature, pore water transformed into cubic ice, which was found to be more stable than cubic ice in the bulk phase, and its stability increased as the pores became smaller. For the smallest pores studied, the cubic ice was stable up to 200 K.

Iatroductioa For a liquid to freeze, an embryo or a cluster of a critical size needs to be formed before homogeneous nucleation can take place. The size of this cluster depends on the temperature and decreases with a decrease in the temperature. When present in pores, not only is the liquid subjected to capiliary forces but also the size of the pore itself may be smaller than the critical cluster size so that the liquid will not freeze at its normal freezing point and a certain amount of supercooling (well in excess of that sometimes required for bulk liquids) will be required for it to freeze. This is well-established through numerous studies on the behavior of liquids and solids in small pores.'q2 Of particular interest has been the study of pore water for which there is reasonable agreement among the various ~ t u d i e s ~ on- ~the depression in freezing point as a function of pore size. However, there is little agreement on whether the heat of melting of pore ice is the same as that of bulk ice. Results indicahg the heat of melting being higher: 10wer,4.~-~ or the same' have been reported. It has also been suggested that the thermodynamic properties of pore ice, in general, should be the same as that of the bulk ice? The melting transition in porous materials is generally not sharp but is spread over a wide temperature range due to the distribution of pore sizes and/or shapes. Furthermore, if the substrate is hydrophilic, a part of the pore water is usually bound to the pore walls and does not undergo a melting transition. The fraction present as bound water is often not known or not determined precisely. All these factors contribute to the uncertainty assoCiatedwith determining the heat of melting. Restricted geometries are known to shift not only the melting transition to lower temperatures but also the solidsolid transition,'' the transition to the superfluid state,",'* and the glass Issued as NRCC No. 34221. 'Present address: Procter and Gamble Pharmaceuticals, P.O. Box 191, Norwich, NY 13815. I Energy Research Laboratories, Energy Mines and Resources Canada, Ottawa, Ontario, Canada KIA OG1.

transition.13J4 However, in the case of water in porous materials, some exceptional results have been reported. In the bulk phase, lowdensity amorphous ice (Ia) is stable up to about 125 K,I5above which it transforms to cubic ice (IC), which then transforms to hexagonal ice (Ih) above 165 K.I6 Dore et al." have reported that on cooling water in 90-A silica pores Ih transforms to ICat 260 K and to a disordered form below 2 13 K and on heating IC remains stable all the way up to the melting point. Cubic ice in the bulk phase is obtained by heating one of the high-pressure forms of ice or amorphous ice.15J6 The transformation from IC to Ih is exothermic and is thus thermodynamically irreversible. In this sense, the direct transformation of Ih to ICas reported by Dore et al. is quite intriguing and warrants further investigation. Hofer et a1.l8 have reported that on uenching water in a hydrogel, with pores in the size range 5-50 , amorphous ice is obtained, which transforms to ICat 210 K and then to Ib at 265 K. Similar results on water in a hydrogel have also been reported by Wilson and Turner,lgand the direct transition of pore liquid water to the amorphous state on quenching was first suggested by BrzhanSzo Thus, in the case of water, the effect of confinement in pores is to raise the phase-transition temperature. Naturally occurring natural gas hydrates are often found in the compacted sediments at the bottom of the Ocean or buried under permafrost. Preliminary estimates2' indicate tbat these hydrates may be the largest reserve of natural gas (mostly CH,). Current interest in these systems is 2-fold: they may serve as a source of energy; they may also contribute to global warming if they become destabilized and release the CH4 Almost all the laboratory studies on gas hydrates have been conducted in the bulk phase. The models dealing with resource assessment and global warming potential of gas hydrates assume that the thermodynamicproperties of gas hydrates in confined spaces are the same as those in the bulk phase. Even if this may k true to a large extent, the stability region of the hydrates will be expected to be shifted to lower temperatures because of lowering of the melting point of ice in pores. It was reported previouslyz3that

w

0022-3654/92/2096-8S94$03.00/0Published 1992 by the American Chemical Society