Monte Carlo Simulations of High-Pressure Phase Equilibria of CO2

ARTICLE pubs.acs.org/JPCB

Monte Carlo Simulations of High-Pressure Phase Equilibria of CO2
H2O Mixtures Yang Liu, Athanassios Z. Panagiotopoulos, and Pablo G. Debenedetti* Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States

bS Supporting Information ABSTRACT:

Histogram-reweighting grand canonical Monte Carlo simulations were used to obtain the phase behavior of CO2
H2O mixtures over a broad temperature and pressure range (50
C e T e 350
C, 0 e P e 1000 bar). We performed a comprehensive test of several existing water (SPC, TIP4P, TIP4P2005, and exponential-6) and carbon dioxide (EPM2, TraPPE, and exponential-6) models using conventional Lorentz
Berthelot combining rules for the unlike-pair parameters. None of the models we studied reproduce adequately experimental data over the entire temperature and pressure range, but critical assessments were made on the range of T and P where particular model pairs perform better. Away from the critical region (T e 250
C), the exponential-6 model combination yields the best predictions for the CO2-rich phase, whereas the TraPPE/TIP4P2005 model combination provides the most accurate coexistence composition and pressure for the H2O-rich phase. Near the critical region (250
C < T e 350
C), the critical points are not accurately estimated by any of the models studied, but the exponential-6 models are able to qualitatively capture the critical loci and the shape of the phase envelopes. Local improvements can be achieved at specific temperatures by introducing modification factors to the Lorentz
Berthelot combining rules, but the modified combining rule is still not able to achieve global improvements over the entire temperature and pressure range. Our work points to the challenge and importance of improving current atomistic models so as to accurately predict the phase behavior of this important binary mixture.

I. INTRODUCTION The CO2
H2O system is one of the most important geological fluids on earth.1 Accurate description of the phase behavior of this binary system is not only a fundamental scientific problem but is also crucial in many industrial and technological fields, such as separation equipment design in the chemical and energy industries,2 enhanced oil recovery, and CO2-based enhanced geothermal systems.3 Moreover, the system has been receiving increased attention due to the rising interest in CO2 capture and geological storage.4 Phase equilibrium of CO2
H2O liquids plays an important role in the assessment of the long-term behavior of carbon dioxide in deep underground reservoirs.5 The CO2
H2O system is known to exhibit gas-gas immiscibility and has a saddle-shaped critical locus that ends in two critical end points.6 Among the numerous experimental studies that describe CO2
H2O mutual solubility (see refs 7 and 8 for a comprehensive review), only a few cover the high-pressure and high-temperature range as well as the corresponding critical r 2011 American Chemical Society

points. Key references include Takenouchi and Kennedy (1964)9 and Todheide and Franck (1963),10 which are frequently cited by later studies. However the quantitative agreement between these two investigations is unsatisfactory9,11,12 and the discrepancies remain unsolved. For example, later experiments by Sterner and Bodnar13 agree with the data of Todheide and Franck, whereas the results of Takenouchi and Kennedy at high temperatures are supported by Blencoe et al.12 Thermodynamic models have been developed to compute the mutual or one-sided solubility of the mixture.8,14 Agreement with experimental data at low temperatures is generally achieved,3,7,15 but even recent models with adjustable parameters still encounter difficulties in reproducing experimental data at high temperatures and high pressures near the critical region.3,5,15 Received: February 15, 2011 Revised: April 4, 2011 Published: April 29, 2011 6629

dx.doi.org/10.1021/jp201520u | J. Phys. Chem. B 2011, 115, 6629–6635

The Journal of Physical Chemistry B

ARTICLE

In this work, we aim to investigate the phase behavior of the system from a molecular simulation perspective. Compared to traditional thermodynamic models, molecular simulation provides a better insight into the microscopic-level properties and their relation to the macroscopic phase behavior. The performance of several existing water and carbon dioxide models in predicting the phase behavior of the CO2
H2O system at low pressures (up to 200 bar) and medium temperature range (75
C ∼ 120
C) has been previously studied using Gibbs ensemble Monte Carlo simulations.16 The same technique was applied later to the CO2
H2O
NaCl system.17 Lisal et al.18 also used NPT Monte Carlo simulations to test the ability of existing models to predict Henry’s constant for carbon dioxide in water. Those studies highlight the importance of selecting accurate models. In the present work, we perform a comprehensive test of several existing water (SPC, 19 TIP4P,20 TIP4P2005,21 and exponential-6 22 ) and carbon dioxide models (EPM2, 23 TraPPE,24 and exponential-625) in predicting the phase equilibrium of the CO2
H2O system over a broader temperature and pressure range (50
C e T e 350
C, 0 e P e 1000 bar) than covered in previous simulation studies. Using histogram-reweighting Grand Canonical Monte Carlo (GCMC) techniques26 and mixed-field finite size scaling analysis,27 we are also able to obtain the critical points at high temperatures. This paper is structured as follows. Methodological details on the potential models and the simulation approach are provided in section II. Results and comparisons to experimental data are presented in section III. Section IV lists the main conclusions from this work and suggests possible directions for future studies.

II. METHODS A. MODEL SYSTEM The SPC, TIP4P, TIP4P2005, and exponential-6 models were used for water, and the EPM2, TraPPE, and exponential6 models were used for carbon dioxide. The SPC model has three sites and uses a Lennard-Jones potential to describe central oxygen
oxygen interactions and fixed-point charges to represent electrostatic contributions, in common with the other models studied. In the TIP4P and TIP4P2005 models, the charge on the oxygen site is moved by a short distance toward the hydrogen atoms to an “M-site”. In the EPM2 and TRAPPE models, both the carbon and oxygen atoms are represented by partially charged Lennard-Jones beads. The total interaction between two molecules i and j with m and n interaction sites is calculated as the sum of Lennard-Jones and Coulomb interactions: 0 2 ! ! 31 ab 12 ab 6 m X n X qai qbj σ σ ij ij @ 4 þ 4εab
ab 5A ð1Þ uij ¼ ij ab ab 4πε r r rij 0 ij ij a¼1 b¼1 where ε0 is the dielectric constant of vacuum, rab ij is the distance between sites a and b of molecule i and j, respectively, qai is the ab charge on site a of molecule i, and εab ij and σij are the LennardJones interaction parameters between sites a and b. In the exponential-6 models, the fixed-point charge representation is preserved but the central Lennard-Jones potential is substituted by exponential-6 interactions:

8 0 2 " #! !6 31 a b ab ab ab m X n > X q q ε r r 6 > , ij i j ij ij m < @ 4 5A, r ab > r ab þ exp Rab
ij 1
ab ij max, ij ab rm, ij rijab uij ¼ a ¼ 1 b ¼ 1 4πε0 rijab 1
6=Rab ij Rij > > : ¥, ab rijab < rmax , ij ab ab ab where εab ij , Rij , rm,ij, and rmax,ij are the exponential-6 parameters. Table 1 summarizes the intermolecular potential parameters for all of the models studied in this work. For interactions between unlike atoms, the Lorentz
Berthelot combining rules were applied, except for the EPM2 potential, for which the characteristic distance σab ij between unlike sites of carbon dioxide molecules was calculated using the geometric mean in accordance with the original work.23 Therefore, the crossinteraction parameters are expressed as: 8 qffiffiffiffiffiffiffiffiffiffi > < σai σbj for a, b ¼ CCO2 , OCO2 for the EPM2 model ab σ ij ¼ 1 a > : ðσ þ σ bj Þ otherwise 2 i ð3Þ

εab ij ¼

qffiffiffiffiffiffiffiffi εai εbj

ð4Þ

Rab ij ¼

qffiffiffiffiffiffiffiffiffiffi Rai Rbj

ð5Þ

1 rmab, ij ¼ ðrma , i þ rmb , j Þ 2

ð6Þ

ab rmax , ij ¼

1 a b ðr þ rmax , jÞ 2 max , i

ð2Þ

ð7Þ

Most of the simulations were carried out in a cubic box of O O , where σCO is the characteristic size of the length L = 8σCO 2 2 oxygen atom of the carbon dioxide molecule. The effects of system size are discussed in section III. Periodic boundary conditions were applied in the x, y, and z directions. The Lennard-Jones and exponential-6 potentials were truncated at rc = 3σO CO2 and longrange corrections were included. The long-range Coulombic interactions between the charged sites were calculated using the Ewald summation technique. The parameter describing the width of the screening charge Gaussian distribution was set to 5/L, and 518 wave vectors were used in the reciprocal space sum. B. Computational Details. Grand Canonical Monte Carlo (GCMC) simulations were performed for all our calculations. Microstates were generated with a combination of displacement, rotation, insertion, and deletion moves. To accept or reject moves, standard Metropolis criteria28 were adopted. In displaceO . In rotation ment moves, the maximum distance was 0.4σCO 2 moves, a randomly selected molecule was rotated by an amount uniformly distributed between
π/5 and þπ/5 about one of the three space-fixed axes chosen at random. The ratio of probabilities of attempted displacement, rotation, insertion and deletion moves was 1:1:5:5. Attempts to insert and delete carbon oxide and water molecules were made with equal probability. The 6630

dx.doi.org/10.1021/jp201520u |J. Phys. Chem. B 2011, 115, 6629–6635

The Journal of Physical Chemistry B

ARTICLE

Table 1. Geometry and Potential Parameters of the Various Models Used H2O SPCa

a

TIP4P b

CO2 TIP4P2005c

exponential-6d

EPM2e

TraPPEf

exponential-6g

H
O
H (
)

109.47

104.52

104.52

109.5

RO
C (Å)

1.149

1.16

RO
H (Å)

1.00

0.9572

0.9572

1.0668

σO (Å)

3.033

3.05

3.029

RO-M (Å)

0.00

0.15

0.1546

0.0

εO/kB (K)

80.507

79.0

83.2 2.753

1.1433

σO (Å)

3.166

3.154

3.1589

3.1947

σC (Å)

2.757

2.80

εO/kB (K)

78.197

78.020

93.2

159.78

εC/kB (K)

28.129

27.0

29.07

qO (e)


0.82


1.04


1.1128


0.7374

qO (e)


0.3256


0.35


0.3233

RO

12

RO(C)

14

rmax-O (Å)

1.75

rmax-O (Å) rmax-C (Å)

0.6923 0.6292

ref 19. b ref 20 c ref 21. d ref 22. e ref 23. f ref 24. g ref 25.

acceptance ratio was between 20 and 50% for displacement and rotation moves, between 0.01 and 0.1% for insertion and deletion moves in liquid-phase simulations, and between 10 and 30% for insertion and deletion moves at vapor-phase conditions. The system always started with an empty box for vapor-phase simulations, while final configurations from previous runs were used as initial configurations for the liquid-phase simulations.29 A total of 10 million equilibrium steps and 50 million production steps were used for vapor-phase simulations while 10
50 million equilibration steps and 100 million production steps were used for liquidphase simulations. For a few simulations at very high densities, where the acceptance ratio was close to 0.01%, we performed several parallel runs at the same condition and combined all histograms until the noise of the reweighted histograms became acceptably small. Near the critical point, the system sampled both vapor and liquid sides and 200 million production steps were used at these conditions. After the system was equilibrated, data on the number of particles of each species, N1, N2, as well as the energy of the system, E, were collected every 300 to 1500 steps, depending on the run’s acceptance ratio, in preparation for generating histograms. A typical run near critical conditions using 200 million production steps took about 24 h on a 3 GHz Xeon processor. To apply the histogram-reweighting techniques26 in binary systems, histograms near the critical point were combined first to construct a free energy landscape covering a broad range of conditions.29 The combined histogram was then reweighted to other conditions, i.e., lower temperatures and broader range of chemical potentials. When the histogram became noisy, new runs at this condition were conducted and the new histograms were added. The phase coexistence curve was calculated by matching the area under the vapor and liquid portion of the density probability distribution. At temperatures far below the critical point, the system could only sample one side of the phase envelope due to large free energy barriers. In this case, we used input conditions which were slightly outside the predicted phase diagram to avoid metastability problems. Histograms from runs on both the liquid and vapor sides were then added to cover the entire range of densities of interest. Once a sufficient number of histograms was collected and combined, the properties of any state within the covered range could be obtained. To calculate the coexistence pressure, we applied the following thermodynamic relation: βP ¼

ln Ξðμ1 , μ2 , V , βÞ V

ð8Þ

where β = 1/kBT, with kB denoting Boltzmann’s constant and T, temperature; P is the pressure; V is the system’s volume; Ξ is the grand canonical partition function; and μi is the chemical potential of species i. ln Ξ is only known up to a constant. To obtain the additive constant, we performed simulations at low densities and calculated the partition function (within the additive constant). At low densities, ideal gas behavior (PV = NkBT) is expected; thus, the value of ln Ξ should scale linearly with N, with unit slope and with the y intercept yielding the additive constant. The mixed-field finite-size scaling approach of Wilding et al.30,31 was applied to mixtures27,32 for determining the critical parameters. For binary systems, this approach defines an ordering operator M, which is a linear combination of the number of particles N1, N2, and total configurational energy U M  N1
sU
qN2

ð9Þ

where s and q are the field mixing parameters. At criticality, the probability distribution of the scaling parameter x, which is defined as x = A(M
Mc), has a universal form. The nonuniversal parameter A and the critical value of the ordering operator Mc are chosen so that zero mean and unit variance of the scaling parameter distribution are obtained. To obtain the critical parameters for the CO2
H2O system, we adjusted the chemical potential, temperature, and field mixing parameters, s and q, so as to achieve the best fit to the three-dimensional universal distribution.

III. RESULTS AND DISCUSSION Five CO2/H2O models were tested in our work: EPM2/ TIP4P2005, EPM2/TIP4P, EPM2/SPC, TraPPE/TIP4P2005, and exponential-6 CO2/ exponential-6 H2O. Simulation results for these models are listed in the Supporting Information. Figures 1 and 2 show the coexisting compositions as a function of pressure at low temperatures (e250
C) for the H2O-rich and CO2-rich phases, respectively. In Figure 1, the experimental data sets9,10,33 are consistent with each other, and the mole fractions calculated from all models have the qualitatively correct monotonic dependence on pressure. Reasonable agreement between experimental data and simulation results is achieved at low pressures (P < 200 bar) for most of the tested models, but overall the simulation results for most models deviate from experiment at higher pressures. We see that the best model for predicting the CO2 solubility in water at 50 and 150
C is the TraPPE/TIP4P2005 model, since it generally agrees with the 6631

dx.doi.org/10.1021/jp201520u |J. Phys. Chem. B 2011, 115, 6629–6635

The Journal of Physical Chemistry B

Figure 1. CO2 solubilities in H2O at low temperatures. Experimental data of ref 9 (9); experimental data of ref 10 ((); experimental data of ref 33 (b); exponential-6 CO2/exponential-6 H2O model (4); EPM2/ TIP4P2005 model (þ); TraPPEþTIP4P2005 model (); EPM2/SPC model (O); EPM2/TIP4P model (3). Solid lines connecting the experimental data are guide to the eye obtained by fitting data to smooth curves.

experimental data across the range of pressures studied. The EPM2/TIP4P and the EPM2/SPC models exhibit similar predictive ability, in that they are both consistent with experiments at T = 150
C, but both overestimate the mole fraction of CO2 at higher temperatures (T > 150
C) and underestimate this quantity at T < 150
C. The solubility of CO2 is severely overestimated by the exponential-6 CO2/ exponential-6 H2O model, and underestimated by the EPM2/TI4P2005 model, across the entire low temperature range. The gas-phase data are more difficult to assess, due to the lack of agreement between experimental data sets from different laboratories. The composition data of the CO2-rich phase reported by Todheide and Franck are generally more consistent with other experimental studies as well as with the results of thermodynamic models between 100 and 260
C in the low pressure (0 to 600 bar) range (see Figure.1 in ref 3). At even lower temperature, e.g., T = 50
C, the H2O mole fraction data of most of the studies are consistent with Wiebe and Gaddy,34 but the data of Todheide and Franck are much higher.7 These discrepancies notwithstanding, the exponential-6 CO2/ exponential-6 H2O model shows the best overall performance in predicting the vapor-phase composition. At 250
C, this model agrees well with the results of Takenouchi and Kennedy; at 150
C, the model is consistent with the data of Todheide and Franck; at 50
C, the predictions are more consistent with the data of Wiebe and Gaddy. The predictions of the TraPPE/TIP4P2005 and EPM2/TIP4P2005 models are consistent with each other, and they both overestimate the mole fraction of CO2. The EPM2/ SPC model and the EPM2/TIP4P model, again, exhibit “crossing” behavior, in that they underestimate the mole fraction of CO2 at higher temperatures (250
C) and overestimate the mole fraction of CO2 at lower temperatures (50
C). In Figure a of Supporting Information, we compare our results for the

ARTICLE

Figure 2. H2O solubilities in CO2 at low temperatures. Experimental data of ref 34 (f); all the other symbols are the same as in Figure 1.

EPM2/SPC model at T = 100
C and low pressures (P e 200 bar) with the simulations of Vorholz et al.16 (shown in filled circles), in which the EPM2/SPC model is also used. The agreement is satisfactory and our results exhibit less scatter than the previous calculations,16 which used the Gibbs Ensemble technique. Figure 3 shows the phase diagram of the CO2
H2O system at high temperatures where critical points appear. Here again, the experimental data are inconsistent with each other. Blencoe et al.12 reported equilibrium data at T = 300
C which agree well with the results of Takenouchi and Kennedy and discussed possible explanations for the discrepancies with respect to the measurements of Todheide and Franck. At these conditions, the simulation results for all the tested model combinations compare poorly with experimental data. The exponential-6 CO2/ exponential-6 H2O model yields the qualitatively correct shape of the phase envelopes at all temperatures, since the model does show critical points at these conditions. But the critical pressures predicted by this model are much lower than the experimental values, and the mole fraction of CO2 in the liquid phase is higher than in experiments, similar to the low-temperature results in Figure 1. The vapor phase predications of this model are more consistent with the data of Takenouchi and Kennedy. The TraPPE/TIP4P2005 and the EPM2/TIP4P2005 models provide better prediction of the liquid branches at 350 and 300
C but not of the gas phase data, and they are both incapable of predicting the intrinsic critical points at 300
C and below. On the other hand, the EPM2/TIP4P and EPM2/SPC models do not exhibit phase coexistence above T = 275
C. At T = 275
C, they both have critical pressures that are much lower than the experimental values, and they both significantly overestimate the mole fraction of H2O in the gas phase. At T = 275
C, we only report the data for the EPM2/TIP4P and EPM2/SPC models because that is approximately the highest temperature where we see the occurrence of critical points with these models. Furthermore, because the TraPPE/TIP4P2005 and EPM2/TIP4P2005 models fail to predict the critical point and hence fail to yield the 6632

dx.doi.org/10.1021/jp201520u |J. Phys. Chem. B 2011, 115, 6629–6635

The Journal of Physical Chemistry B

ARTICLE

Figure 4. Finite size effects on the computed phase diagram of the CO2
H2O system at T = 325
C, as predicted by the exponential-6 models.

Figure 3. Mutual solubilities of CO2
H2O system at high temperatures. Symbols are the same as Figure 1.

qualitatively correct shape of the phase envelope at T = 300
C, we expect them to predict the phase diagram at T = 275
C even more poorly. One of the reasons for most models to yield poor predictions is their generally poor agreement with experimental vapor pressures. The Supporting Information includes a figure showing the vapor pressure and critical point of pure water predicted by the water models studied in this work (Figure b). The critical points of the exponential-6 and TIP4P2005 models are very close to the experimental value, which could explain the fact that the binary phase diagrams using these models better resemble the experimental results at high temperatures(e.g., T = 350
C). Overall, the exponential-6 model demonstrates the best capability in predicting the vapor pressure of water, and consequently it is the best model in predicting the binary phase diagram in the vapor phase (Figure 2) and near critical region (Figure 3). The TIP4P2005 model, on the other hand, deviates from the experimental curve markedly upon decreasing the temperature. Since both the TIP4P and SPC models have much lower critical temperatures than real water, it is not surprising that the binary phase diagram of these models only exhibit phase coexistence starting from a much lower temperature (T = 275
C) than real water. On the other hand, all of the CO2 models yield reliable estimates of the vapor pressure.24,25 Finite size effects on the computed phase transitions were studied: a typical example is shown in Figure 4 for the exponential-6 CO2/ exponential-6 H2O model at T = 325
C. The phase O O , L = 8σCO , diagrams for three system sizes, L = 6σCO 2 2 O L = 10σCO2, are reported and the error bars are calculated by repeating the same simulation three times with different seeds for generating random numbers. The error bars for some of the data are smaller than the symbol size and therefore are not shown. The statistical uncertainties of the data in other figures are expected to be comparable to the uncertainties shown here. Increasing O O to 8σCO slightly increases the CO2 the system size from 6σCO 2 2 mole fraction in the coexisting vapor and liquid phases, but O ) does not change further increasing the system size (to 10σCO 2 the compositions at phase coexistence, which indicates finite size O . effects are negligible below the critical point beyond L = 8σCO 2

Figure 5. Phase diagram of the CO2
H2O system at high temperatures using modified Lorentz
Berthelot combining rules, eqs 10 and 11, with ξ = 0.9, η = 1; (9) Experimental data of ref 9, 275
C; (2) Experimental data of ref 9, 300
C; (@) exponential-6 CO2/exponential6 H2O model, 275
C; (6) exponential-6 CO2/exponential-6 H2O model, 300
C.

The critical points are obtained by the mixed field finite-size scaling analysis for mixtures,27 which is analogous to the pure fluid analysis30,35(see section II), and some residual finite size effects remain for them as shown in the figure. At large enough system sizes, we would expect the predicted critical pressure scale with box length as L-(d
1/ν), where d is the dimension and ν is the correlation length exponent. However, for the limited number of O is too small to system sizes that we studied here, L = 6σCO 2 exhibit such behavior. We also performed some tests of empirical modifications of the combining rules, using the exponential-6 CO2/exponential O . It has been H2O model and a simulation box size L = 6σCO 2 36,37 that further improvement can be shown in previous studies achieved by introducing one or two independent modification factors to eqs 3 and 4, so that qffiffiffiffiffiffiffiffi a b εab ð10Þ ij ¼ ξ εi εj

σ ab ij ¼ 6633

1 ηðσ ai þ σbj Þ 2

ð11Þ

dx.doi.org/10.1021/jp201520u |J. Phys. Chem. B 2011, 115, 6629–6635

The Journal of Physical Chemistry B

Figure 6. Phase diagram of the CO2
H2O system at T = 300
C using modified Lorentz
Berthelot combining rules, eqs 10 and 11, with ξ = 1 and various η values.

We first fixed η = 1 and varied ξ to look for the greatest improvements of predicative ability at high temperatures, e.g., T = 300
C. Figure 5 shows the results for ξ = 0.9 at T = 300 and 275
C, together with the experimental results of Takenouchi and Kennedy. The locations of the critical points are much better predicted at both temperatures, having been shifted to significantly higher pressures. The compositions of both the vapor and liquid phases at T = 300
C are also better predicted. However, the CO2 mole fraction in the liquid phase is underestimated at T = 275
C. The CO2 mole fraction in the vapor phase also increases significantly, becoming inconsistent with the results of Takenouchi and Kennedy. The effects of parameter η are shown in Figure 6, where the results for η = 0.9, 1.0, and 1.1 are compared. Here we see that as η is decreased, the CO2 mole fraction is shifted to higher values for both the vapor and liquid phases. η = 0.9 or 1.1 can each improve one side of the phase diagram, but there does not appear to exist a η value that can improve the prediction for both coexisting phases.

IV. SUMMARY AND CONCLUSIONS In this work, we have comprehensively tested several existing CO2/H2O models (TraPPE/TIP4P2005, EPM2/TIP4P2005, EPM2/TIP4P, EPM2/SPC and exponential-6/exponential-6) in terms of their abilities to predict the mutual solubility of the CO2
H2O system over a broad range of temperature (50
350
C) and pressure (0
1000 bar). We used histogram reweighting grand canonical Monte Carlo simulations and obtained mixture critical points with the mixed-field finite size scaling approach of mixtures. We found that all of the tested models qualitatively capture the temperature and pressure dependence of the solubility, yet none is capable of reproducing the experimental data over the entire temperature and pressure range investigated here. Different models should be used at different conditions. Specifically, the exponential-6/exponential-6 is the best model in predicting the composition of the vapor phase at T e 250
C, whereas the TraPPE/TIP4P2005 model is the most consistent with experimental data on the liquid side. At T > 250
C, only the exponential-6/exponential-6 model can qualitatively predict the presence of the critical points. This model can also capture qualitatively the correct shape of the corresponding phase envelopes at all temperatures.

ARTICLE

We also found that introducing a modification factor to the Lorentz
Berthelot combining rules improves the predictive abilities of the exponential-6/exponential-6 model at specific temperatures, but it is still difficult to achieve global improvements over the entire temperature range by introducing simply one parameter to the Lorentz
Berthelot combining rules. All of the models we tested have simple fixed-charge interactions, effective pair potentials and commonly adopted mixing rules, all of which could limit predictive capability. Further improvement of the quality of predictions using molecular simulations may result from testing the effect of polarizability, manybody interactions, and different combining rules in predicting the phase behavior of pure systems and mixtures. The limitations of existing models notwithstanding, it is important to extend the application of state-of-the-art molecular simulation techniques to the study of the thermodynamic properties and phase behavior of more complex mixtures, such as the CO2
H2O
NaCl system, that are of relevance to carbon mitigation technologies.

’ ASSOCIATED CONTENT

bS

Supporting Information. Complete fluid-phase equilibrium data for the various models studied, vapor pressure of water predicted by the water models studied, and mutual solubilities of CO2
H2O system at low pressures for the EPM2/SPC model. This material is available free of charge via the Internet at http:// pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT P.G.D. gratefully acknowledges the financial support of BP (Carbon Mitigation Initiative at Princeton University). A.Z.P. acknowledges support from the Department of Energy, Office of Basic Energy Science, under Grant DE-SC0002128. ’ REFERENCES (1) Roedder, E. Rev. Mineral. 1984, 12, 1. (2) Rumpf, B.; Nicolaisen, H.; Ocal, C.; Maurer, G. J. Solution Chem. 1994, 23, 431. (3) Spycher, N.; Pruess, K. Transport Porous Media 2010, 82, 173. (4) Socolow, R. H. Sci. Am. 2005, 293, 49. (5) DePaolo, D. J.; Orr, F. M. Phys. Today 2008, 61, 46. (6) Rowlinson, J. S.; Swinton, F. L. Liquids and liquid mixtures, 3rd ed.; Butterworth Scientific: London, 1982. (7) Spycher, N.; Pruess, K.; Ennis-King, J. Geochim. Cosmochim. Acta 2003, 67, 3015. (8) Hu, J. W.; Duan, Z. H.; Zhu, C.; Chou, I. M. Chem. Geol. 2007, 238, 249. (9) Takenouchi, S.; Kennedy, G. C. Am. J. Sci. 1964, 262, 1055. (10) Todheide, K.; Franck, E. U. Ber. Bunsen-Ges. Phys. Chem. 1963, 67, 836. (11) Bakker, R. J.; Diamond, L. W. Geochim. Cosmochim. Acta 2000, 64, 1753. (12) Blencoe, J. G.; Naney, M. T.; Anovitz, L. M. Am. Mineral. 2001, 86, 1100. (13) Sterner, S. M.; Bodnar, R. J. Am. J. Sci. 1991, 291, 1. (14) Ji, Y. H.; Ji, X. Y.; Feng, X.; Liu, C.; Lu, L. H.; Lu, X. H. Chin. J. Chem. Eng. 2007, 15, 439. 6634

dx.doi.org/10.1021/jp201520u |J. Phys. Chem. B 2011, 115, 6629–6635

The Journal of Physical Chemistry B

ARTICLE

(15) Pappa, G. D.; Perakis, C.; Tsimpanogiannis, I. N.; Voutsas, E. C. Fluid Phase Equilib. 2009, 284, 56. (16) Vorholz, J.; Harismiadis, V. I.; Rumpf, B.; Panagiotopoulos, A. Z.; Maurer, G. Fluid Phase Equilib. 2000, 170, 203. (17) Vorholz, J.; Harismiadis, V. I.; Panagiotopoulos, A. Z.; Rumpf, B.; Maurer, G. Fluid Phase Equilib. 2004, 226, 237. (18) Lisal, M.; Smith, W. R.; Aim, K. Fluid Phase Equilib. 2005, 228, 345. (19) Pullman, B. Intermolecular forces: proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry held in Jerusalem, Israel, April 13
16, 1981; D. Reidel; Kluwer: Dordrecht, Holland, 1981. (20) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (21) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123. (22) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1998, 102, 7470. (23) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (24) Potoff, J. J.; Siepmann, J. I. AlChE J. 2001, 47, 1676. (25) Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z. Mol. Phys. 1999, 97, 1073. (26) Ferrenberg, A. M.; Swendsen, R. H. Phys. Rev. Lett. 1989, 63, 1195. (27) Potoff, J. J.; Panagiotopoulos, A. Z. J. Chem. Phys. 1998, 109, 10914. (28) Frenkel, D.; Smit, B. Understanding molecular simulation: from algorithms to applications, 2nd ed.; Academic Press: San Diego, CA, 2002. (29) Panagiotopoulos, A. Z. J. Phys.: Condens. Matter 2000, 12, R25. (30) Wilding, N. B.; Bruce, A. D. J. Phys.: Condens. Matter 1992, 4, 3087. (31) Wilding, N. B. Phys. Rev. E 1995, 52, 602. (32) Panagiotopoulos, A. Z.; Lenart, P. J. J. Phys. Chem. B 2006, 110, 17200. (33) Wiebe, R.; Gaddy, V. L. J. Am. Chem. Soc. 1939, 61, 315. (34) Wiebe, R. Chem. Rev. 1941, 29, 475. (35) Bruce, A. D.; Wilding, N. B. Phys. Rev. Lett. 1992, 68, 193. (36) Stoll, J.; Vrabec, J.; Hasse, H. AIChE J. 2003, 49, 2187. (37) Vrabec, J.; Stoll, J.; Hasse, H. Mol. Simulat. 2005, 31, 215.

6635

dx.doi.org/10.1021/jp201520u |J. Phys. Chem. B 2011, 115, 6629–6635