2900
J. Phys. Chem. 1995,99, 2900-2906
Solubility of Sodium Chloride in Supercritical Water: A Molecular Dynamics Study S. T. Cui and J. G . Harris* Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139 Received: August 9, 1994; In Final Form: October 27, 1994@ We study the solubility of sodium chloride in steam and supercritical water at temperatures from 450 to 550 "C and pressures from 100 to 300 bar as predicted by the simple point charge water model and a simple model of the ions. We calculated the chemical potential of the sodium chloride dimer in supercritical water using a Kirkwood coupling parameter integration in which we computed via molecular simulation the reversable work required to increase the NaC1-water interactions from some reduced value to their full strength. The contribution to the chemical potential of the NaCl particle with reduced interactions we computed using the Widom test particle method. Such a two stage approach is efficient because the strong interactions between NaCl and water makes the Widom insertions of the NaCl inefficient, and a singularity at the point where the NaC1-water interactions become zero makes Kirkwood integrations to the point where the particle vanishes less accurate. We then used the quasiharmonic approximation to predict the chemical potential of the solid salt at the desired temperature and pressure. We find the predicted solubilities agree well with the available experimental measurements. Simple electrostatic continuum pictures do not properly capture the variation of the NaCl solubility with temperature and pressure because they do not take into account the molecular nature of the solvent.
Aqueous electrolyte solutions are important for their industrial applications and their role in biomolecular systems. Recent developments in supercritical water oxidation technology have spurred increased interest in the structure and phase equilibria of these systems at high temperatures and pressures.' Aqueous salt solutions at these conditions are much less well understood than at ambient conditions. Recently researchers have started examining the structure of supercriticalwater and salt Critical to the development of supercritical water oxidation technology is the understanding of phase equilibria in both subcritical brines and supercritical solutions. Several experimental studies have examined the solubility of sodium chloride, one of the most common byproducts of supercritical water oxidati~n.~-IlUsing an ion-water cluster picture, Pitzer and co-workers have developed several empirical models that fit sodium chloride-water phase equilibria at some supercritical conditions.I2-l7 To develop a molecular level understanding of the behavior of supercritical aqueous electrolyte solutions, it is important to identify the features in the intermolecular potentials which are key to determining phase equilibria. Molecular simulations can accomplish this without any approximations beyond the microscopic potential energy surface. In this paper we use such techniques to explore the predictions of the solubility of a salt in an aqueous solution by a simple intermolecular potential model. Thus our calculations provide both a measure of the suitiblity of the model for phase equilibria predictions and an analysis of the effects which determine salt solubility in supercritical water. A few molecular simulation studies of solubility have been carried out recently for relatively simple systems. As the prediction of phase equilibria generally involves calculation of chemical potentials for the two systems in coexistence, the calculation is relatively complex. To simplify the calculation, people have chosen to study simple systems. Swope and AndersenI8 studied the solubility of an inert gas in water, in which case only the chemical potential of the inert gas in water needs to be determined. Shing and ChungIg studied the @Abstractpublished in Advance ACS Absrructs, February 1, 1995.
solubility of a model naphthalene-C02 system, in which all molecules were modeled as spherical particles. In their study, the chemical potential of the pure solid naphthalene was calculated by thermodynamic integration using the experimental vapor pressures and molar volumes. In this paper, we study the solubility of sodium chloride as represented by a simple pair potential model in water as described by the simple point charge (SPC) model at supercritical temperature and pressures. Our previous studies of ion pairing in supercritical water justify assuming that all ions in our model system appear as NaCl dimers.6 Because the concentration of sodium chloride is so low in the saturated solutions we examine, we can assume that ion pairs never interact with each other. We use molecular dynamics simulations to compute the free energies of the solid and solution phases. We then compare the solubility predicted by our simulations with the experimentally measured values and with two different continuum theories.
Phase Coexistence Determination Formalism. The solubility of sodium chloride is the concentration of sodium and chloride ions in the solution at which the chemical potentials of the ions equal their chemical potentials in the solid at the same temperature and pressure. In a previous publication,6 we found that the sodium and chloride ions are strongly paired in supercritical water at the conditions used in most supercritical water oxidation processes. Therefore, we model the sodium chloride solution as a mixture of sodium chloride pairs and water. Mathematically we can write the coexistence condition as 1
&aCI
i
- d i a C l = PNaC1 - PNaC1
(1)
I where pkacl,pNaC1, and pkaclare the chemical potentials of sodium chloride in the solid, in the solution, and as two isolated fixed ions, respectively. The chemical potential of the solid i relative to the isolated ions, pkac, - pNaCI, we will determine from the quasiharmonic approximation (QHA) as described below. The other term can be broken down into
0022-3654/95/2099-290Q~09.00l0 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 9, 1995 2901
The first term on the right hand side is the chemical potential, relative to the isolated ions, of the bound ion pair in an ideal gas composed of ion pairs at a density e, where e is the molar density of water in the supercritical solution. In computing molar densities, each ion pair counts as a single molecule. The second term is the excess chemical potential of the ion pair in the solution, and X N ~ Cis~ the mole fraction of sodium chloride in the solution. We determined p:acl using the molecular simulation techniques described in the next section. Because the sodium and chloride ions are bonded together strongly, it is appropriate to treat them as a harmonic oscillator in computing p:tcl - pkacl and fix the bond length while computing p:acl. p:tc, is the sum of the translation contribution RT ln(@AN,c?), a rotational term, RT In qrot,a vibrational term, RT In qvlb, and the energy of the dimer at the minimum in the potential energy surface, Ubond. The q’s refer to the respective partition functions. The de Broglie wavelength in the partition function is A N ~ = c ~(h2/27cmk~T)1’2, where h is Planck’s constant, k~ the Boltzmann constant, T the temperature, and m the mass of the NaCl molecule. For the rotational partition function the high-temperature limit certainly applies here, and qrot= 87c21k~T/h2, where I is the moment of inertia.20 qvib = exp(-hv/2k~T)/(l - exp(-hv/kBT)) is the partition function for the vibrational degree of freedom, where Y is the fundamental vibrational frequency.20 When the solubility is so low that pyac, is essentially concentration independent in all unsaturated solutions, then once in the infinite dilution limit is known, the concentration of the sodium chloride in the saturated solution can be readily determined by solving eq 1 and 2 for RT In (XNaC1). Molecular Model. Many of the details about the interaction models and the techniques for dealing with the long-range forces we described in a previous paper.6 We review them here briefly. We employ the SPC water model which has point charges of -0.82 at the oxygen atom and 0.41e at the hydrogen atom centers and a Lennard-Jones interaction with a well depth of 0.648 kJ/mol and core diameter of 3.166 A between oxygen centers.21 The ions interact through a Coulombic plus the Huggins-Mayer potential for the sodium and the chloride ions: 22
(3) where the ij indices refer to the ion types and the parameters-Bo , C,, and e;- are given in Table 1. The interactions between the ions and the water molecules are modeled using the potentials introduced by Pettit and Rossky.22 These involve the Coulombic interactions with the SFT charges and Lennard-Jones interactions with parameters listed in Table 1. The one difference between our parameters and those of Pettit and Rossky is that in our model the Coulombic interactions use the SPC charges. Simulation Techniques. We integrated the equations of motion using Andersen’s RATTLE,a velocity Verlet algorithm with bond length constraint^^^ using a time step of 0.0025 ps. Our constant-temperature and -pressure simulations used the Nose- Anderson approach.24 In the calculation of the excess chemical potential of the sodium chloride in water we fixed the distance between the sodium and the chloride ions to that which minimizes the ion pair‘s energy in the vacuum, 2.293 A. The interactions between the ions and the water molecules and between the water
molecules were neglected for distances greater than 9.5 A. For the fluid states we employed the site-site reaction field method introduced by Hummer et al.25 This method replaces the sitesite interaction energy with an effective interaction potential of Vo(r) = q,qj(l/r
+ r2/2R,3 - 3/2R,)
Vo(r) = 0
r > R,
r 5 R, (4)
where R, is the truncation distance. In the reaction field approximation each dipole also contributes a self energy of - I / 21p12/R2to the total potential energy to account for its interaction with the dipole moment it induces on the continuum outside of the sphere of radius R,. In all simulations of the sodium chloride crystal, we employed Ewald sums to treat the long range interactions. Expressions for energy and pressure from Coulombic systems appear in the work of Lantelme et a1.26 Chemical Potential of the Sodium Chloride Molecule in Water. We found a variation of the Kirkwood coupling parameter method most efficient for computing the chemical potential of the sodium-chloride ion pair in water. In the Kirkwood coupling parameter we gradually turn on the interaction between the extra particle and the system by changing the coupling parameter, which we denote by d. That is if the total contribution of the interactions between the ion pair and the rest of the system is denoted as 4 (the dependence upon the spatial coordinates of all of the particles is not shown for convenience), then the difference between the chemical potential of a particle interacting with the rest of the system through potential function d24 and one interacting through a potential function dI4 is
The bracket denotes an average over the configurations of the system in an NPT ensemble, and the subscript d,N+l indicates that the dynamics are for the N 1 system with the interaction between the solute particle and the system given by d4. When d2 is 1 and d1 is 0, this integral gives the excess chemical potential of the inserted molecule. One difficulty that arises is that the integrand becomes singular as d approaches zero. We avoided this problem by evaluating this integral from , I= 0.1 to I = 1.O and using the Widom insertions to determine the excess chemical potential for a particle with interactions given by d 4 when d is 0.1. We integrated eq 5 using the trapezoid rule and integrands at d = 0.1, 0.3, 0.5, 0.7, 0.85, and 1.00. In the Widom test particle method,28a fictitious particle is inserted at various positions of an N-particle system. The excess chemical potential for a constant volume simulation is given by
+
where as in our description of the Kirkwood integration method, 4 is the interaction between the added molecule and the N-molecule system being simulated. In a constant pressure simulation, the expression is modified to take into account the fluctuations of the volume of the system V, so thatI9
Our implementation of the Widom test particle method used the method proposed by Powles et al.,29 in which the three dimensional volume of the system was divided by a certain
Cui and Hanis
2902 J. Phys. Chem., Vol. 99, No. 9, 1995
TABLE 1: Interaction Parameters between the Sodium Ions, Chloride Ions, and WateP 0 Na
H
0.560 1.505
c1 ' €0 and
0.560 1SO5
2.72 3.55
Na 1.31 2.14
4.089E+4 1.210Ef5
Cl
0.317 0.317
1.210E+5 3.360E+5
101.2 674.4
674.4 7045
0.317 0.317
u are the Lennard-Jones well depth and core parameters. B , 40, and C are the parameters of eq 3 in the text.
TABLE 2: Effect of System Size and Cutoff Distance at 250 bar and 800 Ka
e
A V
-80 -100
1
0
2
' '
0.2
' ' '
I
0.4
' ' '
'
0.6
' ' '
0.8
' ' '
' 1
' '
1
N
R, (A)
(U) (M/mol)
AU (kJ/mol)
calculation method
123 123 123 254
9.50 12.66
-234.0 -237.8 -236.5 -236.7
(3.7) (3.4) (4.0) (3.7)
reaction field reaction field Ewald sum reaction field
9.50
O N is the number of water molecules, R, is the cutoff distance, (U) is the average interaction energy between the ion pair and water molecules at the Kirkwood coupling parameter A = 1.0, and AU is the statistical uncertainty in (U).
1.2
h.
Figure 1. Average potential energy of the inserted water molecule (the integrand in eq 5 ) vs the coupling parameter at temperature of 298 K and pressure of 1 bar.
number of grid lines. The test particle was inserted at these grid points. For our system, we put the center of mass of the molecule on the grid points and selected the orientation of the molecule randomly. As a test of our method, we first calculated the excess chemical potential of pure water at the ambient condition T = 298 K, e = 1.0 g/cm3 using 254 water molecules. We used Kirkwood integration to find the difference in the excess chemical potential between 1 = 0.1 and 1.0, while the using the Widom insertions to compute the excess chemical potential at d = 0.1. Since the interaction energy is 10 times smaller than for a full particle, the insertions can converge rapidly. We used 225 ps simulations for each Widom insertion and 475 ps simulations for each Kirkwood integration point for this test and for our studies of the NaCl dimer. The integrand in the Kirkwood expression appears in Figure 1. Integration over the Kirkwood coupling parameter from d = 0.1 to 1.0 gives -35.1 f 1.2 kJ/mol (f95% confidence interval). The excess chemical potential obtained using Widom insertions at d = 0.1 is 11.O f 0.2 Hlmol. To avoid the unrealistic configurations where a hydrogen is inserted on top of an oxygen, we neglected all configurations for which the distance between the oxygen atoms is less than 1.05 A. We predict the excess chemical potential for pure water at this ambient condition to be -24.1 f 1.2 kJ/mol. The excess chemical potential and the excess Helmholtz free energy per particle are related by pex= Aex/N f ( P v )~ @%)id, where 1 and id refer respectively to the liquid and ideal gas state, and v is the molar volume. Using the calculated pressure of 644 bar, we predict Aex/Nto be -22.8 f 1.2, which is in good agreement with refs30and,31which predict excess Helmholtz free energies of -24.3 and -23.9 kJ/mol for the SPC water, respectively. The slight discrepancy with rePo is likely due to statistical error and rePo's using a slightly larger cutoff distance of 10 A and density of 1.017 g/cm3. For the calculation of the excess chemical potential of a sodium chloride ion pair in water, we used 123 water molecules. To ensure that our system size or our treatment of the long range forces do not significantly affect our predictions, we compared the average interaction energies of the NaCl dimer with the water using different potential truncations and system sizes. Table 2 indicates that the dependence of the predicted solubility upon system size and the treatment of the Coulomb interactions
A
-150 -200 -250
0
0.2
0.4
0.6
0.8
1
1.2
h Figure 2. Average potential energy of the inserted sodium chloride molecule (the integrand in eq 5 ) vs the coupling parameter at temperature of 800 K and pressure of 250 bar.
should be insignificant. The methods produce the same energy to within about 1.5%. If this scales linearly with A, its effect is to increase or reduce the logarithm of the solubility by about 1.5%. We compared the excess chemical potential of NaCl in water computed from the Widom insertions with that computed from the Kirkwood integrations. The behavior of the integrand in the Kirkwood integral is illustrated in Figure 2. At 800 K and 250 bar, using the Widom test particle method for the insertion of the NaCl molecule, for a 475 ps run, we obtained a value of -82.0 f 20. kJ/mol. This is within the statistical uncertainty of the -89.4 f 3.2 kJ/mol obtained using Widom method at d = 0.1 and Kirkwood method for 1 between 0.1 and 1.0. Although the results are consistent to within statistical uncertainty, the Widom method suffers from large fluctuations. These large fluctuations occur because for a system of ions in supercritical water, there is a strong clustering of water molecules around the ions. This clustering is stronger in water than in more weakly polar systems because of the strong electrostatic interactions. The Kirkwood coupling parameter method naturally samples those configurations which include water molecules clustering around the ions. The Widom insertion method rarely samples such clusters because the large contributions to the chemical potential occur when one hits a good cluster that has formed without the ions present. This is rather infrequent. At the supercritical water conditions we used, the clustering is primarily due to the Coulombic attractions between the ion sites and the charge sites of the water molecules. At reduced interaction strength, the clustering effect is smaller.
Sodium Chloride in Supercritical Water
J. Phys. Chem., Vol. 99, No. 9, 1995 2903
TABLE 3: Solubility of the Sodium Chloride in Supercritical Water at Various Temperatures and Pressure@ solubility at 723 K (ppmj pressure (bar)
sim
exp
100 150 200 250 300
3.5 8.3 92.0 156.0 494.0
1.5 12.5 63.6
solubility at 773 K (ppm)
Bom approx
sim
exP
5.3 4.3 3.9 2.3
3.6 9.8 43.0 98.0 229.0
0.9 13.7 31.4 101
solubility at 823 K (ppmj
Born approx
sim
exp 0.9
25.3 20.7 17.7 15.8
6.4 13.0 30.0 51.0 162.0
Born approx 107.0 86.7 74.8 66.9
33.8 98.0
a sim represents the values obtained from the simulation, exp represents the experimental values of Armellini and Tester from refs 7 and 8, the Born approx represents the values obtained from the Born approxjmation
Free Energy of the NaCl Crystal. A number of methods have been proposed for the calculation of the free energy of solid Yip et al.36 compared the relative merits for several of them. Most applications of these methods have been to the simple atomic systems, although some have been applied to a molecular liquid such as nitr~gen.~’In this work, we chose the quasiharmonic approximation (QHA) for its simplicity and computational efficiency. For the NaCl crystals we checked the accuracy of the QHA by comparing it with thermodynamic integration from a low temperature. For the NaCl crystal, the Helmholtz free energy, A, in the QHA is given by,
m
3
0
---c Experiment
Born Approx 0
1
~ , , , , , , , , , , , , , , 1 , , , , 1 , , , , , , , , , $
50
100
150
200
250
300
350
Pressure (bar) Figure 3. Solubility of the sodium chloride in supercritical water at temperature 723 K. The filled circles are the simulation result, the diamonds are the experimental result from ref 5, and the open circles are the prediction of the Born approximation.
where D is the 3(N - 1 ) by 3(N - 1) matrix with the elements
D~~~~ = (a2u/ariaar,),
(9)
we chose it to be reasonably close to the experimental one. The system was then run under constant temperature and pressure for at least 62.5 ps, of which the first 12.5 ps were for equilibration. At the end of the run, we computed the average volume and scaled the particle positions to be consistent with this volume. We then employed a conjugate gradient method40 to minimize the energy with the fixed lattice constants. The diagonalization of the dynamic matrix produces the normal frequencies of the system, which in tum determines the free energy through eq 8. When we applied this procedure to solid argon, it produced results in agreement with other researchers’ predictions for argon.41 We tested the reliability of the QHA for the temperature range used in this study by carrying out a thermodynamic integration from 50 K, where the QHA should be extremely accurate, to 800 K . We employed the well-known formula
and the English subscript refers to the particle number while the Greek letters refer to the three Cartesian coordinates of each particle. U is a total interaction potential of the system. The subscript eq refers to the fact that the derivatives are to be evaluated at the equilibrium lattice sites of particles i and j at lattice constants which are the equilibrium average values at the temperature and pressure being studied. UOis the equilibrium static potential energy of the system. A, = (h2/2nmnkBT)”2 is the de Broglie wavelength of the corresponding atom. vn are the vibrational frequencies of the crystal. This first line of the eq 8 is slightly different from the expression for a onecomponent system to take into account the different masses for acsn/ap = Q (10) sodium and chloride. To eliminate the three translational modes, where /3 = l/k& and ( E ) is the average of the total intemal we fix one particle and omit its entries from the matrix. energy of the system at constant volume. The free energy We calculated the cohesive energy and the lattice parameter difference between two temperatures was obtained from the of the sodium chloride crystal at 298 K and 1 bar as a test of integral the accuracy of the model. We obtained a value of -759.3 kJ/mol for the cohesive energy and 2.85 A for the lattice parameter, which are in good agreement with the experimental values of -764.4 kJ/mol and 2.81 A.38 These values are also in excellent agreement with those computed by Fumi and T o ~ i ~ At ~ 250 bar the QHA predicts the integral from 50 to 800 K is 860.1 per ion, while the thermodynamic integration using for the sodium chloride crystal, -758.1 kJImol for the cohesive energy and 2.86 A for the lattice parameter, using the same molecular simulations predicts a difference of 860.2 per ion. Huggins-Mayer potential. Solubility of Sodium Chloride in Supercritical Water. To apply the QHA, we first determined the average molar From the excess chemical potential of sodium chloride in volume at the desired temperature and pressure using a supercritical water and the chemical potential of the sodium molecular dynamics simulation. The simulation was initialized chloride solid, we determined the solubility at various conditions. by putting the sodium and the chloride ions on the sites of two These are listed in Table 3, along with the experimentally interpenetrating fcc lattices, which generate the sodium chloride determined solubilities. Figures 3-5 plot the solubility vs crystal structure. The initial volume of the system can be chosen pressure at three temperatures, 450, 500, and 550 OC. The arbitrarily, but for quick relaxation to the equilibrium volume, solubility decreases with increasing temperature and decreasing
Cui and Harris
2904 J. Phys. Chem., Vol. 99, No. 9, 1995
o
B&n Approx
P v 0
10
l n m n 8 1 s q I~ ' "8 ' I ' " ' I " " J
a
100 r 0
a
v
10
r
+Experiment Born Approx
50
100
150
200
250
300
350
Pressure (bar) Figure 5. Solubility of sodium chloride at 823 K. The symbols are the same as in Figure 3.
Pressure Figure 6. Lines with symbols are -pEac,/kT of the sodium chloride ion pair in water as a function of pressure at 723 K (long dashes, filled circles, 773 K (solid line, squares), and 823 K (short dashes, triangles). The almost horizontal lines without symbols represent (p pR,,,)lkT shifted by a temperature and pressure independent constant. The line types correspond to the same temperatures as in the -p;a,",,,)/kTplot.
pressure. Almost all of the pressure dependence of the solubility can be attributed to the change in the chemical potential of the NaCl ion pair as the pressure is varied. Figure 6 illustrates this with the temperature and pressure dependence of the chemical
potential of the ion pair in the solution and the temperature and pressure dependence of all of the other quantities that get added to the chemical potential to obtain RT lnXNacl-the difference between the chemical potential of the vapor and the chemical potential of the solid (shifted by a constant to fit on the graph). The slight upward trend in (&& - pkac,) is due to the RT In (e)term in the chemical potential. The figures show that the calculated solubility agrees with the experimental data, except at low pressures. [Although the scale is logarithmic and thus the difference between experimental and calculated solub es may appear large when quoted as a percent, because the chemical potentials appear in the exponent, statistical uncertainty alone can be easily over 50%.] The larger discrepancy at low pressures could be due to either the limitations in our model or the lower accuracy of the experimental data at these state points. One of the important features in the solubility data is that it increases rapidly with pressure, reflecting the rapidly increasing solvation power of water to charged particles. The main source in the statistical error comes from the calculation of the excess chemical potential of the sodium chloride in water. The statistical uncertainty (95% confidence limit) in the contributions from the Kirkwood integral ranges from 2.1 kJ/mol at the highest pressure to 3.4 kJ/mole at the lowest pressure. The statistical uncertainties for each value of the integrand were calculated from the standard deviations of 19 block averages, where each block contains the average over 25 ps. The statistical error in the Widom insertion at A = 0.1 has negligible impact as its total contribution is less than 2% of the total excess chemical potential and the statistical uncertainty in it is less than 9% of it. Thus the statistical uncertainties can change the solubilities by a factor of about 1.4-1.8. The static potential energy VOof the sodium chloride crystal and the well depth of the free sodium chloride molecule dominate the contribution to the chemical potential of the solid and the liquid, respectively. From the calculated cohesive energy, the lattice parameter, and other thermophysical properties, we see that the Huggins-Mayer potential with the Coulombic interactions very well describes the crystal solid. However, since the chemical potential is an exponent in the determination of the solubility, some minor changes in the well depth or the shape of the potential can change the solubility significantly. Another factor which could have a significant effect on solubility is the dipole moment of the water molecule. It is well-known that the permanent dipole moment of the SPC water model (2.27) is larger than the experimental permanent dipole moment of water (1.85 D). This has the effect of causing the calculated dielectric constant to be relatively larger that the experimental dielectric constant. Table 4 compares the dielectric constant of SPC water with that of real water and shows that despite the enhanced dipole moment, the SPC model predicts the equation of state of water rather well. If we decrease the dipole moment of the SPC water molecules , we would decrease the interaction energy between the sodium chloride ion pair and the water molecules, thus increasing the chemical potential of the ion pair, which implies a reduction of the solubility. We showed in a previous publication,6 however, that simply scaling the electric charges or the bond lengths of the SPC water would make the equation of state of the SPC much less satisfactory. Thus the other approximations inherent in the SPC model compensate for this deficiency. To compare our results with the predictions of a solventless continuum picture, we used the Born model4*to estimate the
Sodium Chloride in Supercritical Water
J. Phys. Chem., Vol. 99, No. 9, 1995 2905
TABLE 4 temp
(K)
pressure (bar)
723 723 723 723 723 773 773 773 773 773 823 823 823 823 823
100 150 200 250 300 100 150 200 250 300 100 150 200 250 300
exp density" (.g/cm3)
predicted -densityb (g/cm3)
expC
0.03361 0.05420 0.7873 0.1091 0.1484 0.03050
0.0348(16) 0.0571(12) 0.0815(18) 0.1122(16) 0.1464(38) 0.0314(8) 0.0487(14) 0.0689(10) 0.0923(14) 0.1157(12) 0.0281(8) 0.0443(10) 0.0626( 10) 0.0803(10) 0.0997(16)
1.137 1.251 1.418 1.670 2.072 1.114 1.201 1.317 1.472 1.681 1.097 1.167 1.257 1.370 1.511
0.04808 0.06771 0.08990 0.1153 0.02807 0.04363 0.06043 0.07864 0.09844
€
predicted 1.49(4) 1.69(18) 2.04(8) 2.47(12)
An altemative continuum treatment is to bring the ions together in the vapor phase and insert the ion pair into the dielectric continuum with a dielectric constant E . This avoids the difficulties of the unphysical picture which appears when the NaCl bond length becomes small. The chemical potential change upon insertion into the dielectric continuum is given by the expression for the energy of a dipole moment in a spherical
1.38(6) 1.54(10) 1.75(4) 2.00(6) 1.31(4) 1.46(4) 1.60(8) 1.79(4)
Experimental densities were obtained from: Steam Tables in SIUnits Wasserdamjtafeeln; Grigull, U., Straub, J., Schiebener, P., Eds.; Springer-Verlag: Berlin, 1984. Predicted densities are those obtained from our simulations of pure water. Uncertainties are 95% confidence levels. Values of the dielectric constant from a correlation to experimental data were found in: NBSINRC Steam Tables; Haar, L. Gallagher, J. S . , Kell, G. S . , Eds.; Hemisphere Publishing: New York, 1984.
chemical potential of the sodium chloride in water and then used the chemical potential of the solid from the simulation to obtain the solubility of the sodium chloride in supercritical water. The Born model treats the sodium and the chloride ions as charged hard spheres of radius R N and ~ RCI. The change in the Gibbs free energy change due to the transfer of a charged sphere from the vacuum to a continuum with dielectric constant E is
where Ze is the electric charge, R is the radius of the particle, and E is the dielectric constant of the medium, which we obtained from molecular simulations as in our previous work. In our application of the model, we insert the sodium sphere and the chloride sphere separately and then bring them to the desired distance of d = 2.29 8, in the dielectric continuum. The total change of the chemical potential in the procedure is
We define ion radii as the ion-water distance at which the ion-water interaction becomes zero for the appropriate ionwater geometry. In the Na+-water geometry the sodium ion lies on the bisector of the HOH angle with the hydrogens pointing away from the ion. In the Cl--water geometry the C1- lies on the bisector, but the hydrogens point towards the ion. We obtain RNa = 1.88 8,, and Rcl = 2.80 8,. Substituting these values into eq 13, we calculate the Bom estimate for the excess chemical potential of the sodium chloride dimer in water to obtain the solubilities from the Born model, which appear on Figures 3-5. These figures show that the Born model fails to predict the change of the solubility with pressure. In fact it predicts the wrong trend because it underestimates the important ion-ion potential of mean force in solution. When the bond length is so small that l l d becomes larger than ' I ~ ( I / R 4N ~llRCl), the treatment of the ions as hard spheres in a dielectric continuum provides a physically inconsistent picture.
where R is the radius of the cavity. An appropriate radius is given by R3 = CRNaRcl where c = (RNa RCI 412 with d as the distance between the ion centers. This effective radius is 2.637 8,. This approach predicts solubilities which are too small-by a factor of 0.004 for the 723 K, 150 bar state point, and a factor of 0.06 for the 823 K, 300 bar state point. The solubilities increase much too rapidly with pressure using this approximation. Alternatively one can try to fit a radius, R, to match the excess chemical potential of the NaCl in the SPC water using the dielectric constant of the SPC water. Using these properties at the 723 K, 300 bar state point, we obtain a radius of 2.48 A, while using the 150 bar state point at this same temperature we predict a cavity radius of 2.19 8,. The former is close to the mean of the two ion radii, which could serve as another altemative definition of the cavity size. The chemical potential corresponding to the smallest effective radius of 2.11 8, is at the 150 bar, 823 K state point (we did not carry out this test on the 100 bar state points). Overall this effective radius increases with decreasing temperature or increasing chemical potential. Because the cube of the radius is proportional to the cavity volume, the difference between the largest and smallest effective radii (2.48 and 2.11 8,) corresponds to approximately a 60% change in the chemical potential.
+ +
Conclusions We have carried out a study of the solubility of the sodium chloride in supercritical water by molecular simulation from the knowledge of the interaction potential between particles. The calculated solubility is in good agreement with experimental data. The solubility shows a rapid increase with increasing pressure and decreasing temperature. This trend is primarly due to the dependence upon the thermodynamic state of the excess chemical potential of the NaCl ion pair in water. Simple continuum pictures in which the ions are treated as hard spheres in a continuum dielectric do not reproduce the simulation's predictions. The neglect of the strong hydration effect of water molecules around the ion pair is an important reason for the overestimate of the chemical potential in the Born approximation.
Acknowledgment. The author would like to thank Prof. J. W. Tester for useful discussions and D. Lacks for providing to us some of his Monte Carlo simulation results on argon. This work is supported by the Army research Office under the University Research Initiative program. We also acknowledge the Pittsburgh Supercomputer Center for supercomputer time allocations to carry out this work. References and Notes (1) Tester, J. W.; Holgate, H. R.; Armellini, F. J.; Webley, P. A,; Killilea, W. R.; Hong, G . T.; Barner, H. E. In Emerging Technologies in
2906 J. Phys. Chem., Vol. 99, No. 9, 1995 Hazardous Waste Management 111; ACS Symp. Ser. No. 518; American Chemical Society: Atlanta, GA, 1991; Chapter 3. (2) Mountain, R. D. J. Chem. Phys. 1989, 90, 1866. (3) Cummings, P. T.; Cochran, H. D.; Simonson, J. M.; Mesmer, R. E.; Karabomi, S. J. Chem. Phys. 1991, 94, 5606. (4) Cummings, P. T.; Chialvo, A. A,; Cochran, H. D. Chem. Eng. Sci. 1994, 49, 2735. ( 5 ) Gupta, R. B.; Johnston, K. P. lnd. Eng. Chem. Res., submitted. (6) Cui, S. T.; Hams, J. G. Chem. Eng. Sci. 1994, 49, 2749. (7) Armellini, F. J.; Tester, J. W. In 20th lntersociety Conference on Enivornmental Systems; SAE: Williamsburg, VA., 1990; Technical Paper Series 901313. (8) Armellini. F. J.: Tester, J. W. Fluid Phase Eauilib. . 1993, 84, 123. (9) Morey, G. W. Econ. Geo. 1957, 52, 225. (10) Ravich, M. I. Russ. J. Inorg. Chem. 1970, 15, 1041. (11) Sourirajan, S.; Kennedy, G. C. 1962, 260, 115. (12) Bischoff, J. L.; Pitzer, K. S. 1989, 289, 217. (13) Hovey, J. K.; Pitzer, K. S.; Tanger, J. C. J. Phys. Chem. 1990, 94, 1175. (14) Pitzer, K. S.; Pabalan, R. T. 1986, 50, 1445. (15) Tanger, J. C.; Pitzer, K. S. J. Phys. Chem. 1989, 93, 4941. (16) Anderko, A.; Pitzer, K. S. Geochim. Cosmochim. Acta 1993, 57, 4885. (17) Anderko, A.; Pitzer, K. S . Geochim. Cosmochim. Acta 1993, 57, 1657. (18) Swope, W. C.; Andersen, H. C. J. Phys. Chem. 1984, 88, 6548. (19) Shing, K. S.; Chung, S. T. J. Phys. Chem. 1987, 91, 1674. (20) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976. (21) Berendsen, H. J. C.; Postma, J. P. M.; van Gunseteren, W. F. In Zntermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981: Vol. 1; p 331. (22) Pettit, B. M.; Rossky, P. J. J. Chem. Phys. 1986, 84, 5836. (23) Andersen, H. C. J. Comput. Phys. 1983, 52, 24. (24) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384.
Cui and Harris (25) Hummer, G.; Soumpasis, D. M.; Neumann, M. Mol. Phys. 1992, 77, 769. (26) Lantelme, F.; Turq, P.; Quentrec, B.; Lewis, J. W. E. Mol. Phys. 1974, 28, 1537. (27) Straatsma, T. P.; Berendsen, H. J. C.; Postma, J. P. M. J. Chem. Phys. 1986, 85, 6720. (28) Widom, B. J. Chem. Phys. 1963, 39, 2808. (29) Powles, J. G.; Evans, W. A. B.; Quirke, N. Mol. Phys. 1982, 46, 1347. (30) Hermans, J.; Pathiaseril, A.; Andersen, A. J. Am. Chem. SOC. 1988, 110, 5982. (31) Quintana, J.; Haymet, A. D. Chem. Phys. Lett. 1992, 189, 273. (32) Hoover, W. G.; A. C. Hindmarsh; Holian, B. L. J. Chem. Phys. 1972, 57, 1980. (33) Bennett, C. H. J. Comp. Phys. 1976, 22, 245. (34) Frenkel, D. J. Chem. Phys. 1984, 81, 3188. (35) Valleau, J. P.; G. M. Torrie, e. B. J. B. In Modern Theoretical Chemistry; Bome, B. J., Ed.; Plenum: New York, 1976; Vol. 5. (36) Lutsko, J. F.; Wolf, D.; Yip, S. J. Chem. Phys. 1988, 88, 6525. (37) Meijer, E. J.; Frenkel, D. J. Chem. Phys. 1990, 92, 7570. (38) Pauling, L. The Nature of the Chemical Bond and the Structure of molecules and Crystals; an lntroduction to Modern Structural Chemistry; Come11 University Press: Ithaca, New York, 1960; p 526. (39) Fumi, F. G.; Tosi, M. P. J. Phys. Chem. Solids 1964, 25, 31. (40) Press, W. H.; Teukolsky, S. A,; Vettering, W. T.; Flannery, B. P. Numerical Recipes, 2nd ed.; Cambridge University Press: Cambridge, 1992; p 963. (41) Lacks, D. J., personal communications. (42) Born, M. Z. Phys. 1920, 1, 45. (43) Neumann, M. J. Chem. Phys. 1985, 82, 5663. (44) van Gunsteren, W. F.; Berendsen, H. J. C.; Rullmann, J. A. C. Faraday Discuss. 1978, 66, 5 8 .
JP942 109X
