Ion association in concentrated NaCI brines from ambient to supercritical conditions: results from classical molecular dynamics simulations

Highly concentrated NaCl brines are important geothermal fluids; chloride complexation of metals in such brines increases the solubility of minerals and plays a fundamental role in the genesis of hydrothermal ore deposits. There is experimental evidence that the molecular nature of the NaCl–water system changes over the pressure–temperature range of the Earth's crust. A transition of concentrated NaCl–H2O brines to a "hydrous molten salt" at high P and T has been argued to stabilize an aqueous fluid phase in the deep crust. In this work, we have done molecular dynamic simulations using classical potentials to determine the nature of concentrated (0.5–16 m) NaCl–water mixtures under ambient (25°C, 1 bar), hydrothermal (325°C, 1 kbar) and deep crustal (625°C, 15 kbar) conditions. We used the well-established SPCE model for water together with the Smith and Dang Lennard-Jones potentials for the ions (J. Chem. Phys., 1994, 100, 3757). With increasing temperature at 1 kbar, the dielectric constant of water decreases to give extensive ion-association and the formation of polyatomic (NanClm)n-m clusters in addition to simple NaCl ion pairs. Large polyatomic (NanClm)n-m clusters resemble what would be expected in a hydrous NaCl melt in which water and NaCl were completely miscible. Although ion association decreases with pressure, temperatures of 625°C are not enough to overcome pressures of 15 kbar; consequently, there is still enhanced Na–Cl association in brines under deep crustal conditions.


Introduction
The compositions of fluid inclusions 1 indicate that extremely concentrated NaCl brines can be involved in the transport of metals and the formation of hydrothermal ore deposits. There is a variety of evidence that Na and Cl ions become highly associated into NaCl ion pairs in high-temperature aqueous solutions at pressures below 2 kbar. Under such conditions, water and NaCl appear to behave like ideal mixtures 2,3 so that the activity of water, a H 2 O equals its mole fraction (X H 2 O ). Above 4 kbar, however, there is evidence that NaCl ion pairs become dissociated in NaCl-water mixtures. For example, Quist and Marshall 2 find an increase in the electrical conductance of NaCl-H 2 O mixtures between 2 and 4 kbar at 800 uC. More recently, Aranovich and Newton 4 measured the activity of water in concentrated NaCl solutions at 600-900 uC and 2-15 kbar using the brucite-periclase equilibrium. At high pressure, the activity of water is well approximated by a H 2 OX H 2 O /(2 2 X H 2 O ). This is the athermal activity expected if the NaCl-water system behaved as a hydrous salt melt consisting of water and completely dissociated NaCl. 5 The much smaller activity of water under these conditions means that a free solution phase could coexist with an H 2 O-saturated granitic melt in the deep crust. 5 Such a chloride-rich solution phase must play a fundamental role in scavenging metals prior to the formation of hydrothermal ore deposits.
Atomistic simulations of the NaCl-water system could confirm the interpretation of the experimental data and explore PT regimes that are not easily accessible in the laboratory. The simplest atomistic approach is the ''restricted primitive model'' where the solvent is replaced by a dielectric continuum and the ions are approximated as hard spheres with equal radii. Restricted primitive model simulations done by Pitzer and Schrieber 6 and Gillan 7 were used by Oelkers and Helgeson 8,9 to derive association constants of polynuclear clusters in 1 : 1 electrolytes in supercritical aqueous solutions up to 4 kbar. These calculations predicted considerable ion pairing and clustering of Na-Cl in supercritical aqueous fluids. However, the use of a simple constant dielectric to model the solvent means we cannot understand any of the short-range physical interactions. It would be better to use an explicit molecular representation of the solvent. Koneshan and Rasaiah 10 have shown that the degree of ion association in 1 m NaCl (at 683 K with a solvent density of 0.175 g cm 23 ) is signficantly different for continuum and molecular models simulations of the solvent.
A molecular treatment of the solvent and ion interactions is practical for large systems if we can express the interatomic interactions in terms of two-body potential functions. We can use molecular dynamics (MD) or Monte Carlo methods to sample the statistical mechanics of the system and determine not only thermodynamic quantities but also structural information and ionic speciation.
The SPC/E (Extended Simple Point Charge) water model of Berendsen et al. 11 gives a very simple parameterization for the effective charges and short-range potential for water. In the SPC/E model we represent water molecules by three point charges in a rigid geometry with an HOH bond angle slightly different (109.5u) from that in the gas-phase H 2 O molecule (104.5u). Although the model does not allow any polarization or dissociation of the water molecules, it accurately reproduces the thermodynamic and dielectric properties of water, at least along the liquid-vapor coexistence curve. 12 In particular, the SPC/E model predicts the critical point of water to be at 362-374 uC with a density of 0.29-0.326 g cm 23 . For real water, 13 the critical point is at 374 uC with a density of 0.322 g cm 23 .
The SPC/E model also gives a good prediction of the dielectric constant of water and its temperature dependence: the Paper predicted values of e(SPCE)~81.0 at 25 uC and e (SPCE)~6. at T(c) vs. 78.0 and 5.3, respectively, in real water. It is important that the dielectric properties are reproduced as accurately as possible if we wish to predict ion clustering NaCl electrolytes.
Smith and Dang 14 derived ion-water interaction parameters for Na and Cl by that are internally consistent with the SPC/E water model by fitting potential functions to the gas-phase ion-hydration energies. The Smith and Dang 14 parameters (Table 1) are used in a Lennard-Jones potential (for the shortrange interactions) together with a coulombic term: Formal ionic charges q i of 11 and 21 are used for Na 1 and Cl 2 . The parameters were used in molecular dynamical simulations of a NaCl ion pair in 216 water molecules. This was a fairly small simulation but it gave very accurate predictions of experimental hydration energies and complex geometries. Smith and Dang 14 predict that at 25 uC, Na 1 cations will be surrounded by 5.8 water molecules with a Na-O bond length of 2.33 Å . The Cl 2 anions are solvated by 6.9 waters with a Cl-H distance of 2.22 Å . These are in good agreement with those observed experimentally. 15 We are assuming that the interatomic potentials we are using can describe the NaCl-water system at the densities of interest at high pressure and temperature. Physically this means that we are assuming that the bonding and charge distributions in the water molecules at 15 kbar, 625 uC (p~1-1.3 g cm 23 ) are reasonably similar to those at ambient conditions. More sophisticated calculations based on quantum mechanical (ab initio) molecular dynamics are needed to explore this problem. 16 At present, however, such calculations are limited to very small systems for very short run times and cannot address the problem of NaCl ion association as done here. Classical calculations 17 for the equation of state of water at high P and T give good predictions of the water equations of state. 18 This suggests that the potentials used here should be reasonably valid at density conditions of deep crustal fluids.
The formation of ion-pairs and clusters in NaCl solutions has been studied by atomistic MD simulations of dilute to 1 molar NaCl solutions 10,[19][20][21][22] Few simulations, however, have addressed the nature of concentrated NaCl solutions at elevated pressure and temperature. Brodholt 23 obtained pair distribution functions and compared the densities of the simulated solutions to those obtained from the equations of state of Archer 24 and of Anderko and Pitzer 25 and found very good agreement. This was done for solutions up to 5.22 m in concentration and between temperatures of 300 to 2573 K, and under pressures of 1 bar up to 5 kbar. Although Brodholt 23 noted the existence of polynuclear clusters, their existence was not directly tested in his study. The objective of the work presented here is to understand the nature of extremely concentrated (up to 16 m) NaCl solutions at conditions from ambient to those of the lower continental crust (625 uC, 15 kbar).

Computational method and parameters
Molecular dynamics simulations 26 were done on systems consisting of 975-1100 H 2 O molecules and 10-160 NaCl pairs. We believe our systems are large enough to describe the system insofar as simulations on much smaller systems (200-300 H 2 O, 5-20 NaCl pairs) gave essentially the same speciation and density. The simulations were conducted using the MOLDY molecular dynamics program. 27 MOLDY implements periodic boundary conditions and the link-cell method, and uses quaternions to describe molecular orientations. It also uses the Ewald summation in the evaluation of coulombic interactions, and the link-cell algorithm for evaluation of the interactions. The real (r C ) and k-space (k C ) interaction cut-offs and Ewald parameter (a) were determined using eqn. (2)-(4), such that an accuracy of E~exp(2p)~10 25 in the coulombic potential energy was achieved: where N is the number of molecules, V is the volume of the cell and t R /t F has been empirically determined 27 to be 5.5 for MOLDY. For a typical box length of 33 Å (e.g., in a simulation of 1100 H 2 O molecules and 80 NaCl pairs at 1 kbar), this gave r C~1 2.5 Å . We used a time step of 0.25 fs as larger steps caused the system to be unstable. Each run was equilibrated for 25 ps (for temperature scaling, discussed below) and statistics were accumulated for 100-150 ps. Statistics were collected on a constant NPT ensemble. The pressure was kept constant using the method of Parrinello and Rahman, 28 and fixing all the off-diagonal components of the stress matrix to zero. During the equilibration phase, the temperature of the simulation was controlled by re-scaling the velocities every 20 time steps by the factor s~ffi ffiffiffiffiffiffiffiffiffiffi ffi gk B T 2SkT s where g is the number of degrees of freedom, k B is the Boltzmann constant, T the desired temperature, and vkw is the rolling average of the system kinetic energy over the previous 20 time steps. The scaling was carried out separately for translational and rotational components. The re-scaling was turned off once the temperature of the system has reached equilibrium about the desired value, in order to generate a realisation of the constant NPT ensemble. During the unscaled runs, temperatures were constant to within y3%. Fig. 1 shows comparisons of the simulated densities with those obtained from the empirical equations of state. 24,25 There is a systematic and reproducable density discontinuity in the SPC/E simulation that corresponds to the minimum in the excess volume of mixing between NaCl and H 2 O. The NaCl concentration where this occurs is a function of pressure and temperature.

Structure of water in NaCl brines
For pure water, the SPC/E model predicts a local structure at 25 uC, 1 bar (Fig. 2) that results from hydrogen bonding. The resulting O-O pair distribution (Fig. 3) function is in good agreement with results from neutron diffraction. 29 Zhu and Robinson 30 found, in a simulation of a 1.79 m NaCl solution, that the water structure resulting from hydrogen bonding begins to be broken down by the solvation of the Na 1 and Cl 2 ions. At supersaturation (8 m NaCl) we find the local water structure is broken down to force water molecules to occupy the interstitial sites in the hydrogen bonding framework. The change in water structure must explain the saturation concentration of NaCl at these conditions. (Presumably, this simulation would eventually crystallize out NaCl; however, it would take simulation times so long as to be computationally infeasible.) At 325 uC and 1 bar, however, the rotational motion of the water molecules breaks down the long-range structure of pure water; only the O-O peak at 2.8 Å shows up in the pair distribution function of pure water. At the same time, the structure of water in concentrated brines is the same as in pure water. The O-O pair distribution functions do not change with NaCl concentration under PT conditions corresponding to those of the lower crust (Fig. 3b,c).

Na and Cl solvation and ion-association
The coordination numbers and bond lengths for dilute NaCl solutions at 25 uC, 1 bar are in good agreement with those from previous calculations 14 and experiment. 15 For dilute NaCl, we predict a coordination number of 5.3 with only waters in the solvation sphere for Na 1 at 25 uC, 1 bar. With increasing NaCl molality, we find increased Cl complexation of Na. At 25 uC and 1 bar the Na 1 ions in an 8 m solution of NaCl have 4.3 waters and 1.0 Cl ligands in their first coordination shell. There is some evidence of solvent-separated Na-Cl ion pairs from the peak in the Na-Cl pair distribution function at 5 Å . However, this disappears with increasing temperature.
For a given NaCl concentration, the degree of Na-Cl association increases at the expense of Na-O association as we approach the critical point (Fig. 4). For example, the Na 1 ions in an 8 m NaCl solution have, on average, 3.0 H 2 O and 1.8 Cl 2 in their first neighbor shell at 325 uC, 1 kbar. What is especially significant is that there is still considerable association of Na and Cl in concentrated brines under deep-crustal conditions (e.g., at 625 uC and 15 kbar). In the 8 m solution, we have 3.0 waters and 2.0 Cl 2 ions in the first coordination shell. A completely random coordination environment would have 0.6 Cl 2 and 4.4 waters for this composition if the Na 1 coordination number is 5. The ion association persists even at very high concentrations: In the 16 m NaCl solution at 625 uC, 15 kbar, each Na ion is coordinated to, on average, 2.6 Cl 2 ions and 2.7 water molecules. A completely random coordination environment would have 1.2 Cl 2 and 4.1 waters if the total coordination number of Na is 5.3.
The Na and Cl ions form complex polynuclear clusters that are not readily apparent from the pair distribution functions. In order to identify ion pairing and the formation of polyatomic clusters, the coordinates of all the atoms were recorded every 100 timesteps (0.025 ps). These positions were then used to enumerate ion pairs and polyatomic clusters using the assumption that ions within a certain distance R are associated with one another. This distance R was determined from the first minimum in the Na 1 -Cl 2 pair distribution function at 1 kbar, 325 uC and 5 m concentration. The position of the first minimum was roughly constant regardless of temperature or concentration. The resulting speciation is shown in Fig. 5. Previous work 22 has found ion-pairs with short lifetimes in simulations of 1 m solutions at 20 uC. However, we also find   large polynuclear clusters as we approach supersaturation with NaCl. At 325 uC we find somewhat more extensive ion-pairing and the formation of polynuclear clusters (i.e., Na n Cl m with n,m w 3) but the results do not appear qualitatively different from those at 25 uC. Up to NaCl concentration of 3 m, the most important clusters which form are the same as those proposed by Oelkers and Helgeson,8,9,31,32 i.e., NaCl 0 , NaCl 2 2 , Na 2 Cl 1 and Na 2 Cl 2 0 . Other, more complex, clusters (i.e., Na n Cl m with n,m w 3) form with increasing total NaCl concentration. Each of these higher-order clusters has a very small lifetime; collectively, however, they comprise the bulk of the ionic speciation. Species with n,m w 3 exist for less than 5 ps, with some existing for just 0.1 ps or less. In very concentrated brines at high temperature, nearly all of the Na 1 ions are associated with larger clusters. However, there is no significant peak in the Na-Na pair distribution functions that would result from large discrete Na n Cl m clusters resembling nanoscale nuclei of crystalline NaCl. A snapshot of 8 m NaCl at 325 uC and 1 kbar (Fig. 6) shows that the Na-Cl association gives a threedimensional network results form the avoidance of Na-Na and Cl-Cl nearest-neighbor pairs.
A chemically significant consequence of the formation of low-charge or neutral polyatomic clusters is a drop in the ionic strength of NaCl brines at high concentrations. The calculated ionic strengths under ambient conditions and at 325 uC, 1 kbar (Fig. 7) show that concentrated brines have effective ionic strengths less than half of their stoichiometric values.

Dynamics of Na-Cl association
Since conductance measurements have been used to measure ion association in NaCl brines, it is of interest to determine the mobility of ions and clusters as a function of pressure,   temperature and concentration. In the dilute limit, the limiting conductance (L 0 ) of a monovalent electrolyte is where e is the electronic charge, F is the Faraday constant, k B is Boltzmann's constant and D 1,2 are the diffusion constants for the individual ions. The self-diffusion constant (D) of a species can be determined from the Einstein relation where the term in brackets is the average displacement of the ion as a function of time (t). We cannot separate out the mobilities of the different clusters or free ions. However, by determining the concentration effect on the overall selfdiffusion constant for Na or Cl, we can at least see a qualitative effect. The formation of large (Na m Cl n ) n2m clusters should lower the average diffusion constant of Na and Cl if the residence time of an ion in a cluster is long compared to the time between atomic collisions. Calculated self-diffusion constants are shown in Fig. 8. In agreement with that found by Lyubartsev and Laaksonen, 33 the self-diffusion constant decreases with concentration. However, we find that the effect of concentration is relatively small compared with that of temperature and is especially small at the most extreme conditions (625 uC, 15 kbar). Hence, the time-averaged ionic associations are misleading in that they do not show that the clusters have very short lifetimes. Sharygin et al. 34 have recently shown that the high-temperature (378 uC) conductivity of supercritical NaCl solutions at low density need not be fit with models using that include polyatomic NaCl clusters.

Conclusions
The simulations show that with increasing temperature there is increased association of Na 1 and Cl 2 at the expense of solvation by water molecules. This is expected in terms of the decreased dielectric constant of water. The effect of pressure up to 15 kbar is not enough to overcome the effect of temperature to 625 uC. Strong association between Na 1 and Cl 2 ions persist to extreme (625 uC, 15 kbar) conditions. There is little fundamental difference between the structures of the NaCl-water mixtures at 1 kbar, 325 uC and 15 kbar, 625 uC. The molecular picture of fully dissociated ions in high pressure NaCl-water mixtures that has been put forward by previous authors, 4 therefore, is not correct. However, the large clusters have short lifetimes and the residence time of a Na or Cl ion in a polyatomic cluster is comparable to the time between atomic collisions. Because of this dynamic nature of the ion association, the transition from a conventional electrolyte to a hydrous molten salt is continuous with pressure and temperature.