Ion clustering in aqueous salt solutions near the liquid/vapor interface

Molecular dynamics simulations of aqueous NaCl, KCl, NaI, and KI solutions are used to study the effects of salts on the properties of the liquid/vapor interface. The simulations use the models which include both charge transfer and polarization effects. Pairing and the formation of larger ion clusters occurs both in the bulk and surface region, with a decreased tendency to form larger clusters near the interface. An analysis of the roughness of the surface reveals that the chloride salts, which have less tendency to be near the surface, have a roughness that is less than pure water, while the iodide salts, which have a greater surface affinity, have a larger roughness. This suggests that ions away from the surface and ions near the surface affect the interface in opposite ways.


Introduction
Ions in aqueous solutions have different propensities to be near the interface with the vapor phase [1][2][3][4][5][6][7][8]. The interest in the surface affinity follows from the relevance to the environmental properties of the interface and from the connection to interfaces with electrodes and proteins [3,9]. The ions which have a greater affinity for the surface are larger, softer anions such as iodide and thiocyanate, while smaller, harder ions (flouride or alkali cations) tend to avoid the interface. The affinity for the surface is believed to be driven not only by size but by polarizability [1][2][3]. The asymmetric solvation environment of the interface (relative to the isotropic bulk) would favor the formation of a larger induced dipole on the ion, promoting the stability of polarizable ions at the interface. Recent studies have suggested a separate mechanism for surface selectivity, in which the ions stabilize the surface entropically by enhancing fluctuations of the interfacial height [8,10,11]. The ions with a surface affinity have a looser solvation shell which induces interfacial fluctuations larger than either ions with tightly bound solvation shells or the pure liquid.

Methods
The charge transfer model. The simulations use a recent force field which includes both charge transfer and polarizability [29,42,52]. Charge transfer is treated using the discrete charge transfer (DCT) method, in which charge is transferred between pairs based on the distance. For water, a small amount (−0.02e) of charge is transferred from the hydrogen bond acceptor to the donor. For anions, charge is transferred to water based on the distance between the ion center and the hydrogen on the water, and for cations, it is based on the distance between the ion center and the oxygen atom. Charge is transferred between unlike ions based on the distance between ion centers and no charge is transferred between ions of the same charge. The charge transfer amounts are chosen to reproduce the results of quantum chemical calculations. Polarizability is treated in the water molecules using the fluctuating charge formalism, in which charge can redistribute among atoms on the same molecule in response to the electric field due to other atoms, and given a charge constraint determined by the amount of charge transferred with other particles [52]. For ions, polarizability is treated using a charge on a spring, or Drude, model [48]. The models have Lennard-Jones interactions between non-hydrogen atoms and Coulombic interactions between charge sites.
The water-water part of the potential was optimized to reproduce a number of properties of the liquid, including energy, density, pair correlation functions, dielectric constant, and diffusion constant [52]. As relates to this study, the model also accurately reproduces the surface tension [46]. The ion-water potential was parameterized to reproduce single ion solvation free energies, coordination structure, average ion dipole and charge, in the liquid phase, plus ion-water dimer properties [42,48]. Finally, the ion-ion interactions were adjusted to reproduce the osmotic pressure as a function of concentration [29]. The osmotic pressure is sensitive to the amount of ion pairing and provides a good property to optimize ion-ion interactions against [27,28,[53][54][55].
Simulation details. The simulations use 2840 water molecules and 104 ions, creating a 1 M solution.
An orthorhombic box of dimensions 44 × 44 × 176 Å, periodic in all directions. This creates a water layer about 44 Å thick, with a 132 Å thick vapor layer. Ewald sums were used for the electrostatic interactions, with a correction to mimic a system periodic in 2 dimensions [56]. Additional details of the simulations are the same as previously published, [29,48] in the T V N ensemble, at 298 K, using a Nosé-Hoover thermostat and bonds constrained using SHAKE [57]. The four different salt solutions were simulated for a total of 6 nanoseconds. The pure liquid, with 2840 water molecules, in a 44 × 44 × 176 Å box was simulated for comparison to the salt solutions.
Data analysis. The interface was characterized by fitting the water density as a function of the zcoordinant (the direction perpendicular to the interface) to a hyperbolic tangent function, where ρ L and ρ V are the liquid and vapor phase densities, respectively, z 0 is the position of the Gibbs dividing surface, and d is the width. Interfacial widths are commonly reported as 10-90 thicknesses, t l , the length over which the density changes from 90% to 10% of ρ V , which is equal to 2.197d . Properties of the interface are calculated relative to z 0 . Another method to define the surface uses the instantaneous surface (INS), which accounts for the roughness of the surface [58]. This method takes a specific configuration of the system and generates a density using a Gaussian convolution of the heavy atoms. A position of the interface on a three dimensional grid can be found from the location where the density drops to one half of its bulk value. In this analysis, the ions as well as the water molecules are used to define the density [48,59], and a Gaussian width of 2.4 Å and a grid spacing of 1.0 Å is used, as in previous studies [48,58,59]. The INS method allows for characterization of the interface in terms of the fluctuations of the interfacial height [8,10,11]. If the average height of the grid points corresponding to the INS is 〈h〉, the fluctuations in the height can be found from 〈δh 2 〉 = 〈(h − 〈h〉) 2 〉. The shape of the interface can also be characterized by the area of the INS [12]. A useful quantity, A excess , is the area of the instantaneous interface, A, divided by the area of the flat interface, or A excess = A/(L x L y ), where L x and L y are the box lengths in the x and y directions. The ions are grouped into clusters based on pair distances [29]. Ions are taken to be paired if they are within the cut-off distance, r cut , (3.5 Å for NaCl, 4.0 Å for KCl and NaI, and 4.5 Å for KI) and are grouped into larger clusters if any ion is part of more than one pair [34,36].

Results
Properties of the interface. The surface tension of the interface is found from the pressure tensor, as where L z is the box length in the z-direction (perpendicular to the interface) and p αα is the diagonal elements of the pressure tensor in the α direction. Values of the surface tension are given in table 1.
For comparison, the surface tension of pure water using TIP4P-FQ+DCT model [46] is given. The 10-90 thickness, t l , is about the same for pure water, aqueous NaCl and KCl, and increases for KI and NaI. From the instantaneous interface analysis, the fluctuations in the interfacial height, 〈δh 2 〉 1/2 , and the excess surface area, A excess can be determined. The results show that for the NaCl and KCl solutions, 〈δh 2 〉 1/2 and A excess are less than they are for pure water and for the NaI and KI solutions, 〈δh 2 〉 1/2 and A excess are about equal to the pure water values. The 10-90 thickness, t l , is about the same for pure water, aqueous NaCl and KCl, and increases for KI and NaI. This would be consistent with more iodide than chloride ions at the surface. Ion densities. The ion densities relative to the Gibbs dividing surface are shown in figure 1. Iodide shows more surface affinity that chloride, and the cations occupy positions beneath away from the surface. Figure 2 shows the density profiles relative to the instantaneous interface. When viewed this way, the density profiles are sharper and the water density shows some structure [48,59,66]. The anions and cations are seen to be in distinct layers relative to the instantaneous interface.
Ion pairing. The distribution of ion pairs is shown in figure 3. Shown is the distribution of pairs, based on the center of mass of the pair, divided by the total number of the pairs in the whole system. Enhanced pairing at the interface is apparent for all salts, as is especially evident from the distribution relative to the instantaneous surface. Pairing with iodide is enhanced more than chloride, as might be expected since iodide has a greater surface affinity. The increased amount of pairing may arise since there are more anions at the surface or because there are properties of the surface that promote pairing, as suggested by Venkateshwaran, et al. [12].
Ion clustering. To analyze how clustering changes at the surface, the system is split up into three regions, as suggested by figures 2 and 3. The z-coordinate of the center of the cluster relative to the instantaneous surface, z cluster , is used. The first region (z cluster < 5 Å) corresponds to the surface, the second (5 Å < z cluster < 12 Å), the subsurface region, and the third (z cluster > 12 Å) corresponds to the bulk region. Within each region, cluster distributions, p(n), are determined. As done previously, [29] the 23002-3   if all ions are in a n-particle cluster. The distributions are shown in figure 4, comparing the surface to the bulk region. Larger clusters are seen with KI, consistent with the earlier study on bulk solutions [29]. There are differences between the bulk and the surface region, with less larger clusters at the surface.
On the scale of the plot, differences for small clusters are hard to discern, so the p(n) values are given in table 2. KCl and NaI show more, while NaCl shows slightly less and KI significantly less pairing at the surface.    water molecules, while for cations, the charge is relatively unchanged [48]. Pairing decreases the charge of the anion and increases the charge of the cation, as demonstrated with the DCT models [29], ab initio analysis of NaCl structures from classical simulations [51], and ab initio molecular dynamics studies of LiF in water [68]. Pairing leads to a less charge transfer, from the DCT perspective, because there is less charge transfer between an ion pair than there is between an ion and water, particulariy for the anions. At the surface, both an increased pairing and a loss of some of the solvation shell, for anions, leads to a less charge transfer.

23002-5
In the bulk, the charge of the anions is around −0.8e and the cations is 0.9e, leading to an overall charge of the water molecules. This arrises because more charge is transferred to the water from chlo-  ride or iodide than is transferred from sodium or potassium. At a 1 M concentration, there is about 55 water molecules for every ion pair, so the charge imbalance of −0.1e gets among 55 water molecules, giving each an average charge of about −0.002e. This charge on the water molecules, also seen in other simulation studies, [29,42,51] has implications on the dynamics of the water [69]. At the surface, the water molecules can acquire a charge due to hydrogen bond imbalances [44][45][46][47][48][49][50]. Molecules right at the surface tend to accept more hydrogen bonds than they donate to other molecules (often termed "dangling hydrogen bonds") while molecules right beneath the surface donate more than they accept. This gives rise to water molecules which are more negative near the surface (around z = 2 Å) and molecules in the next layer (around z = 6 Å) that are more positive than bulk waters.

Conclusion
The simulation results find that the surface tension depends most strongly on the anion, with a smaller surface tension for iodide. This is in agreement with experiment, which finds that the surface tension increases as KI ≈ NaI < KCl < NaCl [60]. The increases are larger for the model, over 10 dyn/cm for NaCl for example, than they are experimentally (1.64 dyn/cm). The surface tension for pure water varies by over 30 dyn/cm for different water models [46,61,62], so surface tension is very sensitive to the details of the intermolecular interactions. Simulations using non-polarizable models have found good agreement for the surface tension for aqueous NaCl [63][64][65], while the results using polarizable models are not as good, showing a decrease in the surface tension [65]. We are encouraged that the DCT models give the correct trend.
As in previous studies, with DCT models at dilute concentrations [48] and other polarizable potentials [7], the iodide and chloride ions both show an affinity for the surface, with the iodide ion showing more than chloride (figure 1). In dilute solutions, sodium and potassium are repelled from the surface [48]. In the 1 M solutions, the presence of the anion at the surface leads to a peak in the cation distribution beneath the surface, so that there is an interfacial region that is neutral. This peak is slightly enhanced in the KI solution. The ion peaks become more narrow and higher and the water density has the structure, when they are plotted relative to the instantaneous interface (figure 2), consistent with the previous studies [48,59,66], the ion peaks become more narrow and higher and the water density shows some structure. From this analysis, it is apparent that the chloride ion surface peak matches the water surface peak while the iodide peak is shifted towards the vapor phase. The sodium and potassium peaks are shifted more into the liquid phase than they appear in figure 1. In either analysis of the density, there is a region of depleted anion density beneath the surface, consistent with other studies [64,67]. For dilute anions, there is no such region [48], indicating that this feature is induced by the presence of the other ions.
The four solutions all show more pairing at the surface (figure 3). The amount of pairing is a combination of an increased density of anions and the tendency to pair. That is, there might be more pairs because there are more ions or because the surface promotes pairing, as has been proposed by Venkateshwaran et al. [12]. The fraction for each ion to have no other ions in its solvation shell, making it a cluster of size one, is greater for NaCl and KI, and slightly less for KCl and NaI (table 2). Relative to the amount of ions present, there is less pairing for NaCl and KI. For all four ion solutions, there are fewer larger clusters at the surface than in the bulk region ( figure 4). The tendency to form larger clusters is the greatest for KI, as reported earlier for the bulk [29]. For these solutions, there does not seem to be an enhanced tendency to form pairs, or larger clusters, at the surface, possibly reflecting the differences between dilute and concentrated solutions.
The shape of the interface can be characterized by the fluctuations in the interfacial height, 〈δh 2 〉 1/2 , and the instantaneous surface area, as in previous studies [8,[10][11][12]. Here, we define the excess surface area, A excess , as the instantaneous surface area divided by the surface area of the flat interface. Both these properties are measures of the deviations from a flat interface, 〈δh 2 〉 1/2 would be zero and A excess would be one for a perfectly flat interface. They both increase as the surface gets rougher. For NaCl and KCl, the surface is flatter than the pure solution (see table 1). A decrease in the surface area and in the surface height fluctuations is consistent with a large increase in the surface tension of these two salts.
For iodide salts, the surface is rougher than the chloride salts, and 〈δh 2 〉 1/2 and A excess are equal to the

23002-7
values for pure water. Simulation studies on dilute ions suggest that a single ion [8,10,11] or an ion pair [12] enhance surface fluctuations. Combining those results with those of this study suggests that ions in the bulk decrease the fluctuations, and increase the surface tension, and ions at the interface increase fluctuations, and decrease the surface tension.

Acknowledgements
This work was supported by the National Science Foundation under contract number CHE-0611679.