J. Phys. Chem. B 2009, 113, 14277–14290
14277
Monte Carlo Simulations of the Pressure Dependence of the Water-Acid Gas Interfacial Tensions F. Biscay,† A. Ghoufi,‡ V. Lachet,§ and P. Malfreyt*,† Thermodynamique et Interactions Mole´culaires, FRE CNRS 3099, UniVersite´ Blaise Pascal, 63177 Aubie`re Cedex, France, Institut de Physique de Rennes, UMR CNRS 6251, UniVersite´ de Rennes I, 35042 Rennes Cedex, France, and IFP, 1-4 aV. de Bois Pre´au, 92852 Rueil Malmaison Cedex, France ReceiVed: July 22, 2009; ReVised Manuscript ReceiVed: September 3, 2009
We report two-phase Monte Carlo (MC) simulations of the binary water-acid gas mixtures at high temperature and high pressure. Simulations are performed in the NpNAT ensemble in order to reproduce the pressure dependence of the interfacial tensions of the water-CO2 and water-H2S mixtures. The interfacial tension of the binary water-CO2 mixture is determined from 5 to 45 MPa along the isotherm T ) 383 K. Water-H2S interfacial tensions are computed along one supercritical isotherm (T ) 393 K) in a pressure range of 1-15 MPa. The temperature and pressure conditions investigated here by the MC simulations are typical of the geological storage conditions of these acid gases. The coexisting densities and the compositions of the waterrich and acid-gas-rich phases are compared with experiments and with data calculated from Gibbs ensemble Monte Carlo (GEMC) simulations. I. Introduction The risks associated with climate change have been the subject of much debate in recent years. Today, most experts think that these risks are real and directly linked to the increase of the emission of greenhouse gases, especially CO2. Among the various options for mitigating these emissions, large-scale storage of CO2 into underground formations such as deep saline aquifers or depleted oil and gas fields is a promising option.1 There is also a great interest of injecting H2S (together with CO2) resulting from the treatment of natural gases in geological formations. This reinjection of acid gases could be an attractive alternative to the desulfurization of H2S through the Claus process from economic and environmental viewpoints. The success of such CO2 sequestration operations largely depends on the storage safety, which is related to the caprock ability to prevent acid gas leakage. Indeed, acid gases may escape from the formation through different possible pathways, like abandoned wells, faults, or fractures, or through the caprock, which is a low permeability media saturated with water that lies at the top of the reservoir. This risk of leakage through the caprock is largely governed by the fluid-fluid and the fluid-rock interfacial interactions. One of these interactions is the water-gas interfacial tension, that influences the flow process but also controls the capillary-sealing efficiency:
pc ) pwg - pw ) 2γ cos(θ)/R
(1)
where pc is the capillary pressure in the saturated water caprock, γ is the water-gas interfacial tension, R is the largest connected pore throat in the caprock, and θ is the contact angle. Capillary breakthrough occurs when the overpressure, i.e., the difference between the gas pressure pwg in the reservoir and the water * To whom correspondence should
[email protected]. † Universite´ Blaise Pascal. ‡ Universite´ de Rennes I. § IFP.
be
addressed.
E-mail:
pressure pw in the caprock, exceeds the capillary pressure pc. For instance, if R ≈ 0.01 µm and if the water is assumed to be a good wetting fluid (θ ≈ 0), this capillary pressure is about 10 MPa if γ ≈ 50 mN m-1 but is only about 2 MPa if γ ≈ 10 mN m-1. Proper understanding and modeling of these interfacial properties are thus essential in assessing the suitability and the safety of acid gas storage sites. Water-gas interfacial tensions are thus required for the design of such acid gas reinjection projects. Some experimental values of water-CO2 interfacial tension have been reported by several authors since 1957, and a complete literature survey was proposed in 2002 by Hebach et al.2 The reported data indicate however a lack of values in the case of extreme pressure and temperature conditions. Concerning water-H2S systems, very few data are available in the literature, and they mostly correspond to a gaseous H2S phase,3,4 i.e., at pressures (i a)1 b)1 ni nj N-1 N
)
[( ) ( ) ]
∑ ∑ ∑ ∑ 4εab i)1
j>i a)1 b)1
σab riajb
12
-
σab riajb
6
(4)
where riajb is the distance between force center a in molecule i and force center b in molecule j, εab is the energy parameter of the interaction, and σab is the Lennard-Jones core diameter. ni and nj are the number of force centers in the molecules i and j, respectively. N is the total number of water and acid gas molecules. The LJ parameters for the interactions between unlike sites were calculated using the Lorentz-Berthelot combining rules
εab ) (εaaεbb)1/2
1 σab ) (σaa + σbb) 2
(5)
Water-Acid Gas Interfacial Tensions
J. Phys. Chem. B, Vol. 113, No. 43, 2009 14279
The electrostatic interactions were calculated using the Ewald sum method35-38 from
NC
(1) ulrc (zk)
a)1 b)1
∑∑∑ ∑
∑∑
∑∑∑
nj
ni
×
[ ( ) ( )]
j)1
∑ ∑ εab 1/3
∑
1 Q(h)S(h)S(-h) + 2εoV k*0 1 q q erfc(Rriajb)/riajb 8πεo i a j*i ia b jb qiaqib R 1 qia2 erf(Rriaib) 3/2 8πεo i a b*a riaib 4π εo i a
UELEC )
NC
∑ ∑ xi(zk)xj(zk)
8π ) F(zk)2Vs 3 i)1
σ12 ab r9c
-
σ6ab r3c
(11)
(2) ulrc (zk) )
(6)
NC
πF(zk)Vs
NC
∑ ∑ xi(zk) ∫r
∞
c
i)1 j)1
dr
Ns
∫-r d∆z ∑ [Fj(zk+i) r
i)1
Fj(zk+i-1)]rUij,LJ(r) (12) where erfc(x) is the complementary error function and erf(x) is the error function. R is chosen so that only pair interactions in the central cell need to be considered in evaluating the second term in eq 6. The functions S(h) and Q(h) are defined using eqs 7 and 8, respectively
S(h) )
∑ ∑ qia exp(ih · ria) i
(7)
where Fj(zk) and Vs are, respectively, the density of species j and the volume of slab k. ∆z is defined as the difference z - zk. Ns is the number of slabs between z and zk. rc is the cutoff radius fixed to 12 Å for all of the simulations. Uij,LJ(r) is the intermolecular energy with r being the distance between the two centers of mass of molecules i and j.
a
ni
( )
1 h2 Q(h) ) 2 exp - 2 h 4R
Uij,LJ(r) )
(8)
where the reciprocal lattice vector h is defined as h ) 2π(l/ Lxu,m/Lyv,n/Lzw), where u, v, and w are the reciprocal space basis vectors and l, m, and n take values of 0, (1, (2, ... , (∞. The reciprocal space sum is truncated at an ellipsoidal boundary at the vector |hmax|. The convergence factor R is calculated from 2π/Lx. Since the models used here do not involve any electrostatic intramolecular interactions, the last term of eq 6 vanishes. Long-Range Corrections to the Energy for a Multicomponent System. As the geometry of the system shows a heterogeneity along the axis normal to the interface (z-axis), we calculated the LRC to the repulsion-dispersion energy as a function of zk by splitting the cell into slabs of width δz. The total long-range correction energy ULRC was then calculated by summing up all the local contributions of each slab. The ULRC term was then added in the total energy of the system to be used in the Metropolis scheme. The LRCs to the total energy within each kth slab are defined by two parts39 Nz
ULRC )
∑ i)1
Nz
ulrc(zk) )
∑ (ulrc(1)(zk) + ulrc(2)(zk))
(9)
i)1
where Nz is the total number of slabs in the simulation box. In the case of binary mixtures formed by NC components, the total contribution ulrc(zk)40 becomes
ulrc(zk) )
NC
NC
i)1
i)1
∑ Ni(zk)ui,lrc(zk) ) N(zk) ∑ xi(zk)ui,lrc(zk) (10)
where Ni(zk) and N(zk) are the number of molecular species i and the total number of molecules in the slab k, respectively. xi(zk) is the molar fraction of molecular species i in slab k. The (1) (2) (zk) and ulrc (zk) for a binary inhomogeneous expressions of ulrc system are given by
nj
∑∑ a
b
[( ) ( ) ]
4εab
σab r
12
-
σab r
6
(13)
Computational Procedures. The simulation box was an orthorhombic box of dimensions LxLyLz. The sides are of length Lx ) Ly ) 35 Å for the water-H2S system and 32 Å for the water-CO2 system. The details of the geometry of the system are given in Table 2. Periodic boundary conditions were applied in the three directions. The interactions were truncated at 12 Å. The electrostatic long-range interactions are calculated with the full Ewald summation technique with the maximum reciprocal max lattice vectors parallel to the interface given by |hmax x | ) |hy | max ) 8. The increase in |hz | is required to account for the longrange electrostatic interactions accurately in the direction normal to the interface due to a larger dimension box in this direction. |hmax x | was changed from 30 to 70, depending on the pressure used. The number of reciprocal vectors in the z-direction |hzmax| was adjusted during a pre-equilibration period during which Lz fluctuates considerably. For the smallest pressure investigated here, |hzmax| was set to 70 (see Table 2). To reduce the computational time of the simulation, we decreased |hzmax| for the highest pressures and checked the convergence of the electrostatic interactions in the reciprocal space. MC simulations were performed in the NpNAT ensemble. Each cycle consisted of N randomly selected moves with relative probabilities of 0.5 each: translation of the center of mass of a random molecule, rotation of a randomly selected molecule around its center of mass. Simulations were performed at constant normal pressure. This resulted in an additional step in which the volume of the system is changed by varying the longitudinal dimension Lz of the box in each cycle while keeping the interfacial area A constant. Because the longitudinal dimension of the box changes during the course of the simulation, the profiles along the direction normal to the surface are plotted as a function of the reduced z*k defined by zk/Lz. The initial configuration of the water-acid gas interfacial system was prepared from bulk phase configurations of acid gas and water. The cell parameters of the bulk configurations were set to have the same Lx and Ly dimensions. The initial configurations of acid gas and water were built by placing N molecules on the lattice points of a face-centered cubic structure. The nodes of the lattice where molecules were placed as well
14280
J. Phys. Chem. B, Vol. 113, No. 43, 2009
Biscay et al.
TABLE 2: Number of Water and Acid Gas Molecules as a Function of Pressurea pN (MPa)
〈Lz〉
Nwater
Nacid gas
δLz (%)
2000 2000 2000 2000 2000 2000 2000 2000
250 400 650 650 650 650 650 650
Water-CO2 System 5 289.9 7 287.2 10 244.0 15 238.5 20 186.1 30 137.9 40 135.8 45 134.0
0.3 0.4 0.9 1.1 1.8 2.5 2.4 2.6
2500 2500 2500 2500 2500 2500 2500 2500
50 150 250 500 500 750 750 750
Water-H2S System 1.0 277.8 3.0 262.3 5.0 247.4 7.0 295.3 8.5 226.4 10.1 213.9 12.0 168.5 14.6 124.3
0.6 0.7 1.0 0.4 1.3 1.2 2.3 2.7
a The sides are of length Lx ) Ly ) 35 Å for water-H2S and 32 Å for water-CO2 systems. The 〈Lz〉 value corresponds to the average value calculated over the acquisition phase, and δLz represents the maximum variation (%) between the average value 〈Lz〉 and the smallest or largest value of Lz during the acquisition phase.
as the orientation of these molecules were randomly chosen. MC simulations in the NVT ensemble were first performed on these two bulk monophasic fluid configurations. The final configuration was obtained by surrounding the bulk water configuration by two identical bulk CO2 or H2S configurations along the z axis. Some NVT MC step cycles were carried out to equilibrate the interfacial system. A typical MC run in the NpNAT ensemble consisted of 400 000 cycles for equilibration and 600 000 cycles for the production phase. Simulations were performed at 383 K in the pressure range 1-45 MPa for the water-CO2 mixture and at 393 K from 1 to 15 MPa for the water-H2S system. Standard deviations of the ensemble averages were calculating by breaking the production run into block averages. The number of block averages was adjusted in order to allow for the convergence of the surface tension within each block.33 The number of CO2, H2S, and water molecules is given in Table 2. For the lowest pressures (1-8.5 MPa), the number of CO2 and H2S molecules was reduced in order to avoid unreasonable box dimensions and |hmax z | vectors along the z axis. A typical simulation including the equilibration and production phases as well as the calculation of the surface tension takes about 1 month of CPU time over a single processor. It is the price to be paid if we want to model a significant number of molecules and obtain the convergence of the interfacial tension for these binary mixtures. One of the specificities of our MC methodology was the use of the long-range corrections (LRCs) to the configurational energy in the Metropolis scheme. The total long-range correction energy ULRC was updated after each move of molecular position and was added to the energy of the system to be used in the Metropolis scheme. III. Interfacial Tension The most commonly used methods21,25-28,30 for the surface tension calculation are based upon the mechanical route definition and use the tensorial components of the pressure. The first explicit form expresses the components of the pres-
sure tensor as a function of the derivative of the intermolecular potential. This operational expression was given by Kirkwood and Buff21 and is referred to here as the KB expression (γKB). The definition of Irving and Kirkwood26 (γIK) is based upon the notion of the force across a unit area and takes advantage of expressing the local components of the pressure tensor along the direction normal to the surface. A novel method based upon the thermodynamic definition of the surface tension (γTA) has been recently established by Gloor et al.30 and consists of perturbing the cross-sectional area of the system containing the interface. In what follows, we present the different operational expressions of the surface tension with the corresponding expressions of their long-range corrections. Let us consider a system of N molecules with two planar liquid-vapor surfaces lying in the x,y plane where z represents then the direction normal to the surface. Kirkwood-Buff Relation (KB). The molecular surface tension γKB was first introduced by Kirkwood and Buff21 and makes use of the molecular virial expression to give the following relationship:
γKB )
1 2A
〈
N-1
ni
N
〉
nj
rij · riajb - 3zij · ziajb dU(riajb) 2riajb driajb b)1 (14)
∑ ∑ ∑∑ i)1 j)i+1 a)1
where U is the total configurational energy defined in eq 2 and A is the surface area. We used the LRC developed by Blockhuis et al.41 to correct the surface tension of a single component liquid-vapor interface of pure fluids. However, this expression is no longer valid for the liquid-vapor interfaces of binary mixtures. We replace it by the expression developed by Mecke and Winkelmann.42 This LRC expression takes the following form in the case of an AB mixture
Nz
γKB,LRC )
∑
1 γ (z ) 2A i)1 KB,LRC k
(15)
Nz
)
∑
1 (γA (z ) + 2A i)1 KB,LRC k B γKB,LRC (zk)) (16)
where
A γKB,LRC (z*) k ) -24π
(
∫r∞ dr ∫0π dθ[2( 1r )
12
c
-
( 1r ) ] × 6
r2(1 - 3 cos2 θ) sin θ FA(zk + r cos θ) εAB σA σAB /σA 12 π ∞ 2 24π r dr 0 dθ c εA σAB r σAB /σA 6 2 r (1 - 3 cos2 θ) sin θ FB(zk + r cos θ) (17) r
∫
)]
∫
[(
)
This expression can only be used with binary systems and cannot
Water-Acid Gas Interfacial Tensions B γKB,LRC (zk) ) - 24π
(
σAB /σA r
)] 6
σB /σA r
6
c
[(
εAB σA σAB /σA 2 εA σAB r
)
12
-
r2(1 - 3 cos2 θ) sin θ FA(zk + r cos θ) (18) - 24π
( )]
∫r∞ dr ∫0π dθ
J. Phys. Chem. B, Vol. 113, No. 43, 2009 14281
ε σ
[( )
∫r∞ dr ∫0π dθ εAB σBA 2 c
σB /σA r
12
-
The different electrostatic contributions to the pressure tensor are given in Appendix A for completeness. The appropriate LRCs to the normal and tangential components of the pressure tensor in a multicomponent interfacial system have been derived by Guo and Lu40 within the IK definition and are composed of two parts, as expressed in eqs 23 and 24. (1) (2) pN,LRC(zk) ) pN,LRC (zk) + pN,LRC (zk)
r2(1 - 3 cos2 θ) sin θ FB(zk + r cos θ) (19)
NC
be extended to systems with a third component. Irving and Kirkwood Definition (IK). The method of Irving and Kirkwood (IK)26 expresses the surface tension from the local components of the pressure tensor
) -
NC
∑ ∑ xi(zk)xj(zk) ∫r
2π 2 F (zk) 3 i)1 NC
πF(zk)
c
j)1 NC
∑ ∑ xi(zk) ∫r
1 2
∫-LL /2/2 (pN(zk) - pT(zk)) dz z
z
1 A
〈
N-1
N
i)1
j>i
( ) ( )〉
∑ ∑ (rij)R(Fij)β |z1ij| θ
zk - zi zj - zk θ zij zij
(21)
where I is the unit tensor and T is the input temperature. R and β represent x, y, or z directions. θ(x) is the unit step function defined by θ(x) ) 0 when x < 0 and θ(x) ) 1 when x g 0. A is the surface area normal to the z axis. The distance zij between two molecular centers of mass is divided into Ns slabs of thickness δz. Following Irving and Kirkwood, the molecules i and j give a local contribution to the pressure tensor in a given slab if the line joining the centers of mass of molecules i and j crosses, starts, or finishes in the slab. Each slab has 1/Ns of the total contribution from the i-j interaction. The normal component pN(zk) is equal to pzz(zk), whereas the tangential component is given by 1/2(pxx(zk) + pyy(zk)). Fij in eq 21 is the intermolecular force between molecules i and j and is expressed as the sum of all of the site-site forces acting between these two molecules. ni
Fij )
nj
∑ ∑ (fiajb)
a)1 b)1 ni nj
) -
riajb dU(riajb) r driajb b)1 iajb
∑∑
a)1
dr
(22)
drr3
dUij,LJ(r) dr
∫-rr d∆z[Fj(z) -
Fj(zk)]
(20)
where pN(zk) and pT(zk) are the normal and tangential components of the pressure tensor along the normal to the surface, respectively. The method of Irving and Kirkwood26 (IK) is based upon the notion of the force across a unit area. The pressure tensor is then written as the sum of a kinetic term and a potential term resulting from the intermolecular forces. Whereas the first term is well-defined, the potential term is subjected to arbitrariness because there is no unique way to determine which intermolecular forces contribute to the stress across dA. There are many ways of choosing the contour joining two interacting particles. Irving and Kirkwood26 have chosen as a contour the straight line between the two particles. Other choices are possible and result from the lack of uniqueness in the definition of the microscopic stress tensor. The components of the pressure tensor25,27,28 in the Irving and Kirkwood definition are expressed by
pRβ(zk) ) (F(zk))kBT I +
∞
c
i)1 j)1
γIK )
∞
dUij,LJ(r) (∆z)2 (23) dr
As concerns the tangential pressure, the first term is identical to the first term of the normal component, whereas the second term is expressed by NC
NC
∑ ∑ xi(zk) ∫r
π (2) pT,LRC (zk) ) - F(zk) 2 i)1
∞
c
j)1
Fj(zk)]
dr
∫-rr d∆z[Fj(z) -
dUij,LJ(r) 2 [r - (∆z)2] dr
(24)
The first term of pN,LRC(zk) and pT,LRC(zk) is identical to that used in homogeneous molecular simulations except that it requires the use of a local density F(zk) instead of a scalar density F, whereas the second term takes into account the density differences between the slabs. From these LRC expressions of the pressure components, it is then possible to calculate the LRC part relative to the surface tension. These expressions are very close to those established from interfacial systems of a single component.39
γIK,LRC(zk) ) Vs NC π F(zk) 2 A i)1
NC
∑∑ j)1
xi(zk)
∫r
∞
c
dr
∫
r
d∆z -r
Ns
∑ [Fj(zk+i) i)1
dUij,LJ 2 Fj(zk+i-1)] [r - 3(∆z)2] (25) dr The LRC to the surface tension is obtained by summing up all of the contributions to the local values of each bin and dividing the result by 2. Test-Area Method (TA). The test-area method30 (TA) is based upon a thermodynamic route and expresses the surface tension as a change in the free energy43 for an infinitesimal change in the surface area. This infinitesimal change in the area is performed throughout a perturbation process for which the perturbed system (state A + ∆A) is obtained from an infinitesimal change ∆A of the area A of the reference system. The , L(A+∆A) , Lz(A+∆A)) in the perturbed box dimensions (L(A+∆A) x y systems are changed using the following transformations L(A+∆A) x (A+∆A) 1/2 1/2 (A+∆A) ) L(A) ) L(A) ) Lz(A)/(1 + x (1 + ξ) , Ly y (1 + ξ) , Lz ξ),where ξ f 0. The area (A + ∆A) of the perturbed state is (A) (A) (A) thus equal to L(A) x Ly (1 + ξ), and ∆A is equal to Lx Ly ξ. The
14282
J. Phys. Chem. B, Vol. 113, No. 43, 2009
Biscay et al.
operational expression for the calculation of γ within the TA method is
γTA )
k T
B ∑ lim ξf0 ∆A k
〈 (
ln exp -
×
(U(A+∆A)(zk, r′N) - U(A)(zk, rN)) kBT
)〉
(26)
k,A
〈...〉k,A indicates that the average is carried out over the reference state and the k slabs. U(A+∆A)(zk,r′N) and U(A)(zk,rN) are the configurational energies of slab k in the perturbed and reference states. For the calculation of the energy, we adopt the definition of Ladd and Woodcok44 and choose to assign in the slab centered on zk two energy contributions: one contribution due to the energy between the molecules within the slab and a second contribution due to the energy of the molecules within the slab with those outside the slab. The energy of the slab at position zk is defined as N
Uzk )
ni
N
j*i
a
Hk(zi)ULJ(riajb)
k T
k
〈 (
ln exp -
kBT
∑ γLRC(zk)
(30)
k
(28) LJ (zk) can be expressed as After further elaboration, γKBZ
(2),(A) (zk)) (uLRC
)〉
0
Ns
)
0
LJ ELE (zk) + γKBZ (zk)) ∑ (γKBZ
)
× (2),(A+∆A) (uLRC (zk′))
∂A
k
k Ns
As concerns the LRC to the surface tension calculated from the test-area method, the total LRC contribution is expressed as a function of the total long-range correction energy ULRC (see eq 29). The total value of the tail correction of the surface tension within the test-area formalism is given by B ∑ lim ξf0 ∆Aξ
0
∂Uzk
∑ γKBZ(zk)
)
{
∑( ) Ns
)
Ns
b
Hk(zi) ) 1 for zk - δz < zi < zk + δz 2 2 Hk(zi) ) 0 otherwise
Ns
( )
(27)
where Hk(zi) is a top-hat function with functional values of
γTA,LRC )
∂U ∂A
γKBZ )
nj
∑∑∑∑
1 2 i)1
Boltzmann factor. The surface tension value was averaged over the two directions as (γTA,D - γTA,I)/2. However, when the results differ between the two directions, there is no fundamental basis to justify that the errors cancel when averaging the direct and inverse values. One way to check the TA calculation is to compare the average surface tension with other values resulting from different approaches and to check the impact of the value of ξ on the calculated results. We have found that ξ ) 5 × 10-4 represented a reasonable value16,29,32-34,45 for this calculation. Local Expression of the Surface Tension from the Virial Route (KBZ). In a previous work, we established a local version29 of the surface tension based upon the KB expression. The working expression was obtained from the derivative of the potential with respect to the surface and was referred to as the KBZ method in this paper. The reader is directed to Appendix B for further details. The expression of the surface tension within the KBZ approach in the NpNAT statistical ensemble can be written as
(29)
k
The calculation of the surface tension was carried out in the direct (γTA,D) and reverse (γTA,I) directions. The calculation of the direct direction involves an increase of the surface area (A + ∆A), whereas a decrease of the surface area (A - ∆A) is performed in the reverse path. The ensemble average was carried out over the configurations of the reference system, whereas the configurations of the perturbed system did not participate to the Markov chain of states. Thermodynamic consistency requires that the surface tension in the direct and reverse directions must be equal in magnitude and in opposite sign. This is satisfied when the configuration space of the perturbed system matches the one of the reference system. This requirement implies using an appropriate value of ξ. The value of ξ must satisfy two constraints:30 this value should be small enough to allow an accurate calculation of the surface tension from eq 26 and large enough to provide reasonable statistics for the
LJ γKBZ (zk) ) -
48ε
1 ij A(riajb, rij) ∑ ∑ ∑ ∑ riajb riajb i∈k
a
j*i
b
×
[( ) ( ) ] σij riajb
12
-
1 σij 2 riajb
6
(31)
where A(riajb,rij) is given in Appendix B. The analytical ELE (zk) is also given in Appendix B. The local expression of γKBZ versions γTA(zk) and γKBZ(zk) can be compared to the standard local version γIK(zk) and allow to check the features of the profile of the surface tension calculated from the new approach TA and the well-known approach of Kirkwood and Buff. In fact, the local version of the surface tension using the KBZ approach can be considered as the local version of the Kirkwood and Buff definition. The values of surface tension calculated from the local expressions are identical to those calculated from the macroscopic expressions (of γKBZ and γTA) within a maximum deviation of 0.01%. The working local LRC expression for the surface tension within the NpNAT ensemble using the KBZ definition is given by the following equation
Vs NC γKBZ,LRC(zk) ) πF(zk) 2A i)1
NC
∑ ∑ xi(zk) ∫r
(
×
j)1
∞
c
∫-rr d∆z[Fj(z) - Fj(zk)]
)(
(32)
dr
∂Uij,LJ(r) r2 - 3(∆z)2 Uij,LJ(r) + r r ∂r
)
(33)
Water-Acid Gas Interfacial Tensions
J. Phys. Chem. B, Vol. 113, No. 43, 2009 14283 TABLE 3: Interfacial Tension Values (mN m-1) of Water-CO2 (383 K) and Water-H2S (393 K) Binary Mixtures Calculated from MC Simulations Using Different Operational Expressionsa γKB pN
γIK
γLRC
γ
γLRC
5 7 10 15 20 30 40 45
0.81 0.71 0.61 0.71 0.51 0.41 0.51 0.31
41.824 37.023 33.119 30.220 27.917 24.818 24.720 22.318
1.21 1.11 0.71 0.91 0.61 0.51 0.61 0.51
0.95 3 5 7 8.5 10.1 12 14.6
1.51 0.91 0.91 0.81 0.71 0.61 0.41 0.31
43.729 37.322 32.019 24.118 19.812 15.68 10.36 10.05
1.71 1.41 1.21 1.01 1.01 0.91 0.71 0.51
γTH γ
γLRC
γKBZ γ
γLRC
γ
Water-CO2 43.928 0.81 38.227 0.61 34.225 0.71 31.826 0.91 29.024 0.51 26.122 0.41 25.223 0.51 23.319 0.31
43.024 37.126 33.325 29.722 27.519 24.820 22.920 21.618
0.91 0.71 0.81 0.91 1.01 0.61 0.51 0.41
43.523 37.625 33.425 30.124 28.022 25.721 23.320 22.819
Water-H2S 45.127 1.21 39.019 0.91 33.318 0.81 25.015 0.61 20.312 0.71 16.49 0.51 11.16 0.31 11.26 0.41
43.726 37.627 31.016 23.012 18.314 14.18 9.25 8.95
1.31 1.01 0.91 0.71 0.61 0.51 0.41 0.31
44.420 38.316 32.117 24.118 19.29 15.58 10.05 9.95
a
The long-range corrections (γLRC) and the total surface tension (γ) are given for each method. The subscripts give the accuracy of the last decimal(s); i.e., 41.824 means 41.8 ( 2.4.
Figure 1. (a) Water-CO2 and (b) water-H2S interfacial tensions as a function of pressure. The water-CO2 interfacial tensions are calculated along the isotherm T ) 383 K in the range 5-45 MPa, and the water-H2S interfacial tensions are simulated at 393 K from 1 to 15 MPa. We report the experimental water-CO2 interfacial tension of Chiquet et al.5 and of Shah et al.6 The experimental water-H2S interfacial tensions were measured by Shah et al.7 The water-acid gas interfacial tensions46 predicted from the Cahn gradient theory are also given for comparison.
IV. Results and Discussion Interfacial Tensions. Figure 1a reports the comparison between the calculated and experimental5,6 water-CO2 interfacial tensions at 383 K from 5 to 45 MPa. We also add the interfacial tension values46 calculated from the Cahn gradient theory.47 The interfacial tensions calculated from the different definitions KB, KBZ, IK, and TA are given in Table 3 along with the corresponding LRC values. First, we observe in Figure 1a an excellent agreement between the simulated and experimental interfacial tensions. All of the definitions give consistent values within the statistical fluctuations. The average deviation between the simulated and experimental data5 in the pressure range 5-45 MPa data are 2.4, 3.0, 1.1, and 2.6% with IK, KB, KBZ, and TA, respectively. The average deviation is increased up to 11% when we compare the interfacial tensions calculated with MC with those calculated with the Cahn gradient approach. This figure shows that the interfacial tension decreases significantly with increasing pressure from 5 to 20 MPa. Over this pressure range, this corresponds to an interfacial tension decrease of 36% for MC and 35% for the experimental data.5 From 20 to 45 MPa, γ decreases more slowly for both the calculated and experimental interfacial tension. We also observe that the values calculated from the Cahn gradient theory deviate much more from experiments in the high pressure region, whereas the values predicted from MC simulations follow the experi-
mental trend. Table 1 shows that the LRC to the water-CO2 interfacial tension contributes between 1.5 and 3% to the total value. Figure 1b displays the variation of the water-H2S interfacial tension with pressure from 1 to 15 MPa at 393 K. The calculated interfacial tensions are reported in Table 3 with their LRC contributions within the IK, KB, KBZ, and TA approaches. When the interfacial tension is averaged over the different definitions, we note that the maximum deviation between the average calculated surface tensions and experimental data7 is less than 10%. The average deviation between MC and experimental data is close to 6%. Figure 1b confirms this good prediction and shows that the calculated interfacial tensions match very well with experiments and those6,46 predicted from the Cahn gradient theory. Figure 1b shows that the calculated and experimental water-H2S interfacial tension decreases linearly with pressure from 1 to 12 MPa. From 12 MPa, the interfacial tension reaches a plateau of 10 mN m-1. The MC simulations of the water-acid gas binary mixtures in the NpNAT statistical ensemble give excellent predictions of the pressure dependence of the interfacial tension at typical geological formation temperatures and pressures. Table 4 compares the calculated and experimental interfacial tensions at 10 and 15 MPa between the water-CH4, water-CO2, and water-H2S mixtures. The water-CH4 interfacial tensions were calculated in a recent work20 and averaged over the IK, TA, KB, and KBZ definitions. The experimental γexp mixtures are averaged over the data of Ren et al.,48 of Jennings and Newman,49 of Wiegand and Franck,50 and of Tian et al.51 The results in Table 4 show that the water-H2S interfacial tensions at 10 MPa are about 3.3 and 2.0 times smaller than those of the water-CH4 and water-CO2 mixtures for the same temperature conditions. At 15 MPa, the ratios γwater-CH4/γwater-H2S and γwater-CO2/γwater-H2S are increased by 5 and 3, respectively. It means that the capillary pressure of eq 1 is significantly reduced when H2S is present in the reservoir. As a consequence, the risk of leakage is increased when the injected gas is rich in
14284
J. Phys. Chem. B, Vol. 113, No. 43, 2009
Biscay et al.
TABLE 4: Interfacial Tensions of Water-CH4,20 Water-CO2, and Water-H2S Mixtures at High Temperature and Pressure Conditionsa
Figure 2c confirms the good agreement between the target pressure and the calculated total normal component in the CO2 phase. The average value of the normal pressure calculated over the slabs of the water and CO2 phases is 5.7 ( 1.7 MPa. When the normal pressure is averaged only over the slabs of the CO2 phases, pN is equal to 5.02 ( 0.02 MPa. The analysis of the pressure profiles is required in order to check that the total normal pressure within the system is in line with the target pressure. This is an essential step toward the prediction of the pressure dependence of the interfacial tension. Part d of Figure 2 displays the profiles of the LRC parts to the normal pressure. We see that these profiles correspond to relatively small contributions to the total pressure. However, the calculation of these LRC components to the pressure is relevant because it is also in that of the LRC to the surface tension. Solubilities and Coexisting Densities. Figure 3a shows the molecular density profiles of water and CO2 in the mixture at 393 K and 10 MPa. We observe that the water liquid phase extends over a region of approximately 0.4 reduced units. This corresponds to a real thickness of 100 Å. The bulk gas phases of CO2 shown in Figure 3b are also well developed except at 45 MPa where we only observe a small portion of homogeneous density. Close to the interfaces, an enhancement of density at the water surface is observed. The ratio of the maximum CO2 density at the water interface to bulk gas density is greater at low pressure than at high pressure. The maximum CO2 density is twice higher than the bulk CO2 density at 10 MPa and only 1.2 times higher at 30 MPa. For the water-H2S system, the ratio between the maximum H2S density and the bulk gas density (see Figure 3d) decreases from 4.0 to 1.2 as the pressure increases from 7 to 12 MPa. Parts a and c of Figure 3 show that the water density profiles can be fitted by a hyperbolic tangent function. We have calculated the number of CO2 and H2S molecules per unit area involved in the film formed at the interfacial region. This number is calculated from the integral of the gas density profiles between the two z positions delimiting the adsorption layer divided by the interfacial area. Figure 4 shows this number for the water-CO2 (left axis) and water-H2S (right axis) systems. We observe that this number increases from 0.004 molecules Å-2 for the smallest pressures up to 0.12 molecules Å-2 for the highest pressures investigated here. The evolution of this factor with pressure is identical for the two systems, although the range of pressure is different. These ratios correspond to decreases of 45% for the water-CO2 interfacial tension and of 75% for the water-H2S interfacial tension within the range of pressure investigated here. A correlation between the width of the adsorption layer and the decrease of the interfacial tension can be deduced. For comparison, this number does not exceed 0.033 molecules Å-2 in the case of the water-CH4 mixture20 and corresponds to a decrease of only 20% of its interfacial tension from 1 to 50 MPa. Parts a and b of Figure 5 show the hydration number of water molecules within the first hydration shell of water as a function of z*. k This first coordination shell is characterized by a radius of 3.5 Å. The coordination number of CO2 in the first coordination shell of CO2 characterized by a radius of 6.5 Å is also represented for comparison. We check that the wetting layer is defined by an increase of the coordination number of CO2 molecules in the interfacial region, whereas the coordination number of the water molecules decreases steeply in the interfacial zone. The coordination number of CO2 molecules is dependent on the pressure, whereas the first hydration shell of water is not affected by the pressure.
system
p (MPa)
T (K)
〈γ〉 (mN m-1)
γexp (mN m-1)
∆E (kJ mol-1)
water water-CH4 water-CO2 water-H2S water-CH4 water-CO2 water-H2S
10 10 10 15 15 15
388 373 383 393 373 383 393
56.140 50.825 33.525 15.49 48.420 30.423 10.06
55.9 52.315 34.4 15.3 50.216 30.8 10.5
44.95 41.86 33.86 26.63
a The interfacial tension 〈γ〉 values are averaged over the IK, TA, KB, and KBZ approaches. The experimental γexp interfacial tensions of the water-CH4 mixture are averaged over the data of Ren et al.,48 Jennings and Newman,49 Wiegand and Franck,50 and Tian et al.51 The water-CO2 and water-H2S experimental interfacial tensions correspond to the data of Chiquet at al.5 and of Shah et al.,7 respectively. We add for comparison the surface tension of the liquid-vapor interface of water at 388 K calculated from the average of the different definitions16 with the experimental value.52 The ∆E contribution corresponds to the difference of the total energy between the acid gas or vapor phases and the water bulk phase.
H2S. The MC simulations reported here accurately reproduce the different ratios between the water-acid gas and water-CH4 interfacial tensions. We give for comparison the surface tension of the liquid-vapor equilibrium of water52 at 388 K. We also report in Table 4 the difference in the total energy (∆E) between the vapor or acid gas phases and the water phase. This energy contribution calculated from eq 27 is averaged over the slabs of the acid gas phases and water phases. The profiles of the total energy in each phase are not shown here. This ∆E contribution corresponds in part to the energy for generating the interface. As expected from the values of interfacial tensions, we observe that this contribution decreases from 44.9 kJ mol -1 for the liquid-vapor interface of water to 26.6 kJ mol -1 for the water-H2S system. The decrease of this contribution is in line with that of the interfacial tension. This indicates that the decrease of the interfacial tension with respect to the system studied can be explained by different magnitudes of energy contributions between water and methane, water and CO2, and water and H2S. From a methodological viewpoint, we report in Figure 2a the difference between the normal and tangential pressure profiles at 20 MPa. We report the profiles along the normal to the interface as a function of the reduced quantity z*k ) zk/Lz. This pN(z*) k - pT(z*) k profile considers the sum of the LennardJones and the electrostatic contributions. The LRC contributions are not taken into account in this profile. The electrostatic contributions to the pressure are calculated using eqs 36, 37, and 38. The normal and tangential components are the same in the bulk water and CO2 phases with larger fluctuations in the water phase. The two peaks at the interfacial regions are symmetric, and the contribution from both is the same. The integral of this profile γ(z*) k shown in the right axis of Figure 2a is approximately constant in the bulk water and CO2 phases, indicating that the bulk phases do not contribute to the interfacial tension as expected for a system at mechanical equilibrium with well developed bulk phases. Since the simulations are performed in the NpNAT ensemble, we check in Figure 2b that the sum of the profiles of the different pressure contributions matches with the imposed pressure (pN ) 5 MPa). As expected, we observe larger fluctuations of the total pressure in the water phase, whereas the agreement is excellent in the CO2 bulk phases.
Water-Acid Gas Interfacial Tensions
J. Phys. Chem. B, Vol. 113, No. 43, 2009 14285
Figure 2. (a) pN(z*) k - pT(z*) k profile (20 MPa) in the water-CO2 system for the sum of the Lennard-Jones and electrostatic parts as a function of z*k with its integral ∫0z*k (pN(z) - pT(z)) dz represented on the right axis. (b) Profiles of the ideal, LRC, and configurational parts of the pressure. The short dashed line represents the value of the target pressure (5 MPa). (c) Zoom of the total pressure profile in the CO2 bulk phase. (d) Profiles of the LRC to the pressure.
The water phase can be regarded as nearly incompressible in the pressure range. The same features can be observed on the profiles of the coordination number of water and H2S molecules in parts c and d of Figure 5. Typical configurations of the water-acid gas interfacial systems are shown in the different parts of Figure 6. Parts a and b of Figure 7 show the phase diagram of the water-CO2 binary system at 383 K and the bulk water-rich and CO2 -rich phase densities as a function of pressure, respectively. The experimental solubilities of CO2 in the water-rich phase and the experimental composition of the CO2-rich phase correspond to the data of Takenouchi and Kennedy,53 whereas the experimental densities of the water-CO2 system reported in Figure 7b have been measured by Chiquet et al.5 We check in Figure 7a that the direct twophase MC simulations of the water-CO2 mixture accurately describe the compositions of the water- and CO2-rich phases. The maximum deviation observed with experiments is 12% for the solubility of CO2 in water. The increase of solubility of water in the CO2 phase is 1.5 times higher at 40 MPa than at 10 MPa. The MC simulations predict an increase of 1.6 for this pressure range. We report the solubilities calculated from NpT-Gibbs ensemble Monte Carlo (GEMC)54,55 simulations by Vorholz et al.56 using the TIP4P water model57 and the EPM2 CO218 model. Our MC calculations are in line with their solubilities calculated at low pressures. The coexisting densities of the water- and CO2-rich phases from our MC simulations are compared in Figure 7b with the experimental densities5 and the GEMC calculated densities.56 Our MC results predict very well the bulk water-rich density
with a maximum deviation of 3% from experiments in the 5-45 MPa pressure range. The bulk CO2-rich density is reproduced with an average deviation of 6%. We observe that the water densities calculated from GEMC are understimated with respect to both experiments and our calculations, whereas the CO2 bulk densities (GEMC) match very well with our values at low pressure. The molar fractions of the CO2 in the water-rich and CO2-rich phases are given for each simulated pressure in Table 5 for the GEMC and MC simulations as well as for experiments. Figure 7c shows the phase diagram of the water-H2S binary mixture at 393 K. The experimental solubilities of H2S in water and the composition of the vapor H2S-rich phase58 are reported with those calculated in this work and GEMC calculations.59 The phase diagram is well reproduced at low pressures. At higher pressures, the calculated solubilities are in line with those calculated from GEMC. The GEMC simulations59 used the TIP4P model for water and the Kristof and Liszi potential.19 The deviation between our values and experiments is less than 10% in the low pressure range. The composition of the H2S-rich phase, i.e., the solubility of water in H2S, is accurately predicted at low pressure. At higher pressures, the MC solubilities are line with those calculated from GEMC.59 Figure 7d displays the pressure dependence of the bulk water-rich phase and H2S-rich phase densities calculated from MC and GEMC59 with the experimental data.7 The calculated (MC and GEMC) and experimental molar fractions of the H2S in the water-rich and H2S-rich phases are listed for each simulated pressure in Table 5. The average deviation between our results and experiments is about 3%
14286
J. Phys. Chem. B, Vol. 113, No. 43, 2009
Biscay et al.
Figure 3. Molecular density profiles of (a) water and CO2 for the water-CO2 mixture at 383 K and 10 MPa and (b) only CO2 in the pressure range 10-45 MPa. Molecular density profiles of (c) water and H2S for the water-H2S mixture at 393 K and 7 MPa and (b) only H2S from 7 to 12 MPa. The solid line in parts a and c represents the fitted tangent hyperbolic function of the profile.
than that of the GEMC simulation. This is the price to be paid to obtain both the interfacial tension and the compositions of the two phases in physical contact. V. Conclusions
Figure 4. Pressure as a function of the number of acid gas molecules per unit area (molecules Å2) involved in the wetting layer close to the interfacial region. The ratios in the water-CO2 and water-H2S binary mixtures are represented by open circles and full triangles, respectively.
for the densities of the H2S phase and less than 1% for the water bulk phase. The water densities calculated from GEMC using the TIP4P model are underestimated, as it was already mentioned for the water-CO2 mixture. These results demonstrate that the two-phase simulation methods can reproduce quantitatively the phase diagram of binary mixtures as well as the GEMC simulations. However, the computational cost of a two-phase simulation is greater
Monte Carlo simulations have been performed to calculate the water-acid gas interfacial tensions at high temperature and high pressure. We have simulated the water-CO2 mixture in the 5-45 MPa pressure range along the isotherm T ) 383 K and the binary water-H2S system along the supercritical isotherm T ) 393 K in a pressure range up to 15 MPa. These temperature and pressure conditions for these two acid gases are those encountered in the geological storage. The water-acid gas interfacial tensions have been calculated with the mechanical and thermodynamic definitions. All of the calculated values show that the MC simulations can reproduce the pressure dependence of the interfacial tension of the binary water-CO2 with a maximum deviation of 3% from experiments. The calculation of the binary water-H2S interfacial tension deviates by less than 10% from experiments in the 1-15 MPa pressure range. These simulations demonstrate that the two-phase molecular simulation is an attractive tool for producing accurate interfacial tensions of complex mixtures at high temperature and high pressure. Additionally, the compositions of the water-rich phase and the CO2-rich phase are satisfactorily predicted with a maximum deviation of 12% for the CO2 solubility in water.
Water-Acid Gas Interfacial Tensions
J. Phys. Chem. B, Vol. 113, No. 43, 2009 14287
Figure 5. Profiles of the coordination number along the direction normal to the surface for (a and b) the water-CO2 and (c and d) the water-H2S binary mixtures at different pressures.
calculated densities deviate by less than 3% from experiments. We have added for comparison the compositions and the densities calculated from GEMC by other authors. The calculations performed here cover a wider range of pressure, but the agreement is good for low pressures. Concerning the binary water-H2S mixture, our MC simulations predict very well the compositions of each bulk phase within 10% from experiments, whereas the coexiting densities agree with experiments with an average deviation of 3%. The coexisting densities and the solubilities calculated from GEMC compare very well with our values. The agreement is good even though the GEMC calculations with the TIP4P model underestimate the water phase density. Acknowledgment. The authors would like to acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant SUSHI (Grant No. ANR-07-BLAN-0268) “SimUlation de Systemes He´te´rogenes et d’Interfaces”. A. Expression of the Local Components of the Pressure Tensor Calculated from the IK Definition
Figure 6. Configurations of the (a and b) water-CO2 and (c and d) water-H2S binary mixtures with zooms in the interfacial region.
The coexisting densities of the water and CO2 bulk phases are also well reproduced with respect to experiments. The
In the case of dispersion-repulsion and electrostatic interactions calculated from the Lennard-Jones potential and Ewald summation method, respectively, the components of the pressure tensor become within the Irving-Kirkwood definition LJ R K,1 pRβ(zk) ) (F(zk))kBT I + pRβ (zk) + pRβ (zk) + pRβ (zk) + K,2 (zk) (34) pRβ
14288
J. Phys. Chem. B, Vol. 113, No. 43, 2009
Biscay et al.
Figure 7. (a and c) (P,x) phase diagrams of the water-CO2 (T ) 393 K) and water-H2S (T ) 383 K) binary mixtures. The solid lines correspond to the experimental solubilities of the water-CO253 and water-H2S58 mixtures. The solubilities calculated from MC and GEMC56,59 are represented by circles and triangles, respectively. (b and d) Bulk water and acid gas phase densities as a function of pressure calculated from the MC method (circle) and GEMC (triangle). The experimental densities are given for the water-CO25 and water-H2S7 binary mixtures and represented in dotted lines. The full symbols represents the properties in the water-rich phase, and the open symbols refer to the properties in the acid-gas-rich phase. The GEMC calculations use the TIP4P water model,57 the EPM2 CO2 model,18 and the potential18 proposed by Kristof and Liszi for H2S. The MC calculations use the TIP4P/2005 water model17 and the same potentials as GEMC for CO2 and H2S.
where the superscripts LJ, R, and K represent the contributions of the Lennard-Jones interactions, the Ewald real space, and the Ewald reciprocal space, respectively. The Lennard-Jones contribution to the pressure tensor can be calculated from the following expression: LJ pRβ (zk) ) -
N
ni
(rij)R · (riajb)β duLJ(riajb) 1 θ × riajb driajb |zij | b)1
( ) ( )) zk - zi zj - zk θ zij zij
K,2 PRβ (Zk)
(35)
(
∑
(
〈∑ ∑ ∑ ∑ N-1
N
ni
nj
i)1 j)i+1 a)1 b)1
)
(rij)R · (riajb)β 2 × Rriajb exp(-R2riajb2) + erfc(Rriajb) riajb3 √π
( ) ( )〉
)〉
(37)
〉
(38)
na
∑ ∑ (ria - ri)βq a)1
ia
×
∑ Hk(zi)Q(h)ihR[S(h) exp(-ih · ria) -
S(-h) exp(-ih · ria)]
where Hk(zi) is a top-hat function defined in eq 28 and na the number of atoms in molecule a.
qiaqjb ×
zk - zi zj - zk 1 θ θ |zij | zij zij
〈
N
1 2π )4πε0 V2 i)1 h*0
The contribution of the real space to the local pressure is given by eq 36
1 ) 4πε0A
〈
1 2π H (z )Q(h)S(h)S(-h) × 4πε0 V2 h*0 k i 2hRhβ hRhβ δRβ 2 h 2R2
nj
i)1 j)i+1 a)1
R pRβ (zk)
K,1 pRβ (zk) )
1 × A
(∑ ∑ ∑ ∑ N-1
The two contributions of the reciprocal space are given by eqs 37 and 38
(36)
B. Expressions of the Surface Tension within the NpNAT Ensemble Partition Function. In the case of binary systems, we can consider an ensemble characterized by a constant value of the interfacial area A and a constant normal pressure pN (input
Water-Acid Gas Interfacial Tensions
J. Phys. Chem. B, Vol. 113, No. 43, 2009 14289
TABLE 5: p-x-y Values of the Water-CO2 (383 K) and Water-H2S (393 K) Binary Mixtures Where x and y Denote the Molar Fraction of the Acid Gas in the Water-Rich and Acid-Gas-Rich Phases, Respectivelya MC p 5 7 10 15 20 30 40 45
GEMC
x
y
0.006262 0.008383 0.012390 0.015290 0.018690 0.025890 0.026991 0.027395
0.98020 0.97520 0.96618 0.96319 0.95917 0.95815 0.94116 0.93520
0.007979 0.015891 0.023190 0.040295 0.045998 0.058898 0.060494 0.063799
0.83642 0.91546 0.94447 0.95347 0.95247 0.95346 0.95743 0.92446
y
Water-CO2 0.005915 0.009912 0.012626 0.014130
0.9385 0.9543 0.9653 0.9673
x
y
0.0140
0.986
0.065997 0.062383 0.064198
0.9146 0.9339 0.9516
0.0167 0.0256 0.0359
0.9178 0.963 0.946
(A-1)
(
∫ ∫ drN dV exp - U(r k) B+T pN N
V
))
(A-2)
where Vo is an arbitrary unit of volume that makes the partition function dimensionless, rN is the potential energy, and pN is the normal pressure of the system. Λ is the de Broglie thermal wavelength. The surface tension γ is defined as the thermodynamic derivative of the Gibbs free energy with respect to the interacial area A as
( ∂G ∂A )
γ)
((
NpNAT
N
) -kBT
( ) ∂V ∂Lz
and
∂L ∂V ( ∂A ) · ( ∂V ) · ( ∂L∂A ) z
0.9556 0.9566 0.9158
G ) -kBT ln QNpNAT
γ)
) Lz
Lz
pressure). Lz is the longitudinal dimension of the system and fluctuates in the course of the simulation. The resulting statistical ensemble is NpNAT. In a classical system constituted by N identical particles of mass m defined by their coordinates rN and momenta pN, the Gibbs free energy G is related to the partition function QNpNAT by the following expression:
(
(A-5)
)A
(A-6)
A
The different partial derivatives relating the variables V, A, and Lz obey the following equation:
The GEMC calculations56,59 used the TIP4P,57 EPM2,18 and Kristof and Liszi19 potentials for water, CO2, and H2S molecules, respectively. The experimental mole fractions of the water-CO253 and water-H2S58 mixtures are given for completeness. We only report in this table the GEMC and experimental molar fractions determined at the simulated pressure. All of the other mole fractions are represented in Figure 7. The subscripts give the accuracy of the last decimal(s); i.e., 0.9385 means 0.938 ( 0.005.
1 3N Λ N!Vo
∂V ( ∂A )
0.979 0.976 0.974
a
) -kBT ln
z A
with
Lz
0.0210 0.0240 0.0260
0.017323 0.026143 0.046197
z
Lz
Exp.
x
Water-H2S 0.95 3.0 5.0 7.0 8.5 10.1 12.0 14.6
∂V ( ∂A ) dA + ( ∂L∂V ) dL
dV ) LzdA + AdLz )
1 ∂Q Q ∂A
( )
) ( ))
∂U(r ) ∂V + pN ∂A ∂A
NpNAT
NpNAT
A
) -1
(A-7)
z V
From eq A-7, it turns out that the derivative (∂V/∂A)A does not exist. It then follows that (∂V/∂A)A ) 0. The surface tension of eq A-4 becomes
〈( ) 〉 〈( ) 〉 ∂U(rN) ∂A
γNpNAT ) )
∂U(rN) ∂A
Lz
(A-8)
A
+ pNLz
(A-9)
where the angular brackets denote an average over the NpNAT ensemble. KBZ Formalism. In a previous work, we established a local version29 of the surface tension based upon the KB expression. The working expression was obtained from the derivative of the potential with the respect to the surface and was referred to as the KBZ method in this paper. Within the NpNAT statistical ensemble, the derivative of the Lz dimension with respect to the interfacial area becomes
( ) ( ∂Lz ∂A
)
Lz
∂(V/A) ∂A
)
) Lz
Lz -V + )0 2 A A
(A-10)
Using the result of eq A-10, the expression of the A(u,v) parameter becomes
(A-3)
(
A(u, v) ) (u)x
(A-4)
The calculation of the second term of the right-hand side of eq A-4 requires expressing the derivative of the volume with respect to A according to
(v)x 4Lx
2
+ (u)y
(v)y 4Ly2
)
(A-11)
and the definitive expression of the electrostatic contribution to the KBZ definition in the NpNAT ensemble is given by eq A-12
14290
J. Phys. Chem. B, Vol. 113, No. 43, 2009
ELE γKBZ (zk) )
∑
[
∑∑ i∈k
a
2a
√πriajb
1 Q(h) Im Vε0 h*0
)
{
1 4πε0
∑ ∑ -q q
ia jb
j*i
b
exp(-(R2riajb2)) +
{( ∑ ∑ [
-(ri)x
i
a
A(h, ria) (exp(-ihria))(
[
1 A(r , r ) × riajb iajb ij
erfc(Rriajb) riajb2
]
+
πl πm - (ri)y 2(Lx)3 2(Ly)3
∑ ∑ exp(ihr ))} +
]
+
ia
i
a
] ( )[
1 2π2l2 2π2m2 1 h2 + exp 2V h*0 (L )4 (Ly)4 h4 4R2 x
∑
Biscay et al.
]
h2 + 1 S(h)S(-h) 4R2
}
(A-12)
where Im denotes the imaginary part of the complex variable. References and Notes (1) IPPC Special Report on Carbon dioxide Capture and Storage; Cambridge University Press: Cambridge, U.K., 2005. (2) Hebach, A.; Oberhof, A.; Dahmen, N.; Ko¨gel, A.; Ederer, H.; Dinjus, E. J. Chem. Eng. Data 2002, 47, 1540. (3) Herrick, C. S.; Gaines, G. L., Jr. J. Phys. Chem. 1973, 77, 2703. (4) Strathdee, G. G.; Given, R. M. J. Phys. Chem. 1976, 80, 1714. (5) Chiquet, P.; Daridon, J. L.; Broseta, D.; Thibeau, S. Energy ConVers. Manage. 2007, 48, 736. (6) Shah, V. Ph.D. thesis, Universite de Pau et des pays de l’Adour, 2008. (7) Shah, V.; Broseta, D.; Mouronval, G.; Montel, F. Int. J. Greenhouse Gas Control 2008, 2, 594. (8) MacLeod, D. Trans. Faraday Soc. 1923, 19, 38. (9) Weinaug, C.; Katz, D. Ind. Eng. Chem. 1943, 35, 239. (10) Guggenheim, E. A. J. Chem. Phys. 1945, 13, 253. (11) Zuo, Y.-X.; Stenby, E. H. Can. J. Chem. Eng. 1997, 75, 1130. (12) Rowlinson, J. S. J. Stat. Phys. 1979, 20, 197. (13) Cahn, J. W.; Hilliard, J. E. J. Chem. Phys. 1958, 28, 258. (14) Bongiorno, V.; Scriven, L. E.; Davis, H. T. J. Colloid Interface Sci. 1976, 57, 462. (15) Yang, A. J. M.; Fleming, P. D. I.; Gibbs, J. H. J. Chem. Phys. 1976, 64, 3732. (16) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2008, 128, 154716. (17) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (18) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (19) Kristof, T.; Liszi, J. J. Phys. Chem. B 1997, 101, 5480. (20) Biscay, F.; Ghoufi, A.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2009, 131, 124707. (21) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1949, 17, 338. (22) Buff, F. B. Z. Elektrochem. 1952, 56, 311. (23) MacLellan, A. G. Proc. R. Soc. London, Ser. A 1952, 213, 274. (24) McLellan, A. G. Proc. R. Soc. London, Ser. A 1953, 217, 92. (25) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon Press: Oxford, U.K., 1982.
(26) Irving, J. H.; Kirkwood, J. G. J. Chem. Phys. 1950, 18, 817. (27) Walton, J. P. R. B.; Tildesley, D. J.; Rowlinson, J. S.; Henderson, J. R. Mol. Phys. 1983, 48, 1357. (28) Walton, J. P. R. B.; Gubbins, K. E. Mol. Phys. 1985, 55, 679. (29) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. Phys. ReV. E 2008, 77, 031601. (30) Gloor, G. J.; Jackson, G.; Blas, F. J.; de Miguel, E. J. Chem. Phys. 2005, 123, 134703. (31) Ibergay, C.; Ghoufi, A.; Goujon, F.; Ungerer, P.; Boutin, A.; Rousseau, B.; Malfreyt, P. Phys. ReV. E 2007, 75, 051602. (32) Biscay, F.; Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Phys. Chem. B 2008, 112, 13885. (33) Biscay, F.; Ghoufi, A.; Lachet, V.; Malfreyt, P. Phys. Chem. Chem. Phys. 2009, 11, 6132. (34) Biscay, F.; Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2009, 130, 184710/1. (35) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press: Oxford, U.K., 1989. (36) Smith, E. R. Proc. R. Soc. London, Ser. A 1981, 375, 475. (37) Heyes, D. M. Phys. ReV. B 1994, 49, 755. (38) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. J. Chem. Phys. 1995, 102, 4574. (39) Guo, M.; Lu, B. C. Y. J. Chem. Phys. 1997, 106, 3688. (40) Guo, M.; Lu, C. Y. J. Chem. Phys. 1998, 109, 1134. (41) Blokhuis, E. M.; Bedeaux, D.; Holcomb, C. D.; Zollweg, J. A. Mol. Phys. 1995, 85, 665. (42) Mecke, M.; Winkelmann, J.; Fischer, J. J. Chem. Phys. 1999, 110, 1188. (43) Ghoufi, A.; Malfreyt, P. Mol. Phys. 2006, 104, 2929. (44) Ladd, A. J. C.; Woodcok, L. V. Mol. Phys. 1978, 36, 611. (45) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2008, 128, 154718. (46) Shah, V.; Broseta, D. Langmuir 2007, 23, 12598. (47) Cahn, J. W. J. Chem. Phys. 1977, 66, 3667. (48) Ren, Q. Y.; Chen, G. J.; Guo, T. M. J. Chem. Eng. Data 2000, 45, 610. (49) Jennings, H. Y.; Newman, G. H. J. Pet. Sci. Eng. 1971, 171. (50) Wiegand, G.; Franck, E. U. Ber. Bunsenges. Phys. Chem. 1994, 98, 809. (51) Tian, Y.; Xiao, X.; Zhu, H.; Dong, X.; Ren, X.; Zhang, F. Wuli Huaxue Xuebao 1997, 13, 89. (52) Data taken from the saturation properties of water at the http:// webbook.nist.gov. (53) Takenouchi, S.; Kennedy, G. C. Am. J. Sci. 1964, 262, 1055. (54) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (55) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527. (56) Vorholz, J.; Harismiadis, V. I.; Rumpf, B.; Panagiotopoulos, A. Z.; Maurer, G. Fluid Phase Equilib. 2000, 170, 230. (57) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (58) Kuranov, G.; Rumpf, B.; Smirnova, N. A.; Maurer, G. Ind. Eng. Chem. Res. 1996, 35, 1959. (59) Vorholz, J.; Rumpf, B.; Maurer, G. Phys. Chem. Chem. Phys. 2002, 4, 4449.
JP906953A
