Surface Propensities of the Self-Ions of Water.

The surface charge of water, which is important in a wide range of chemical, biological, material, and environmental contexts, has been a subject of lengthy and heated debate. Recently, it has been shown that the highly efficient LEWIS force field, in which semiclassical, independently mobile valence electron pairs capture the amphiproticity, polarizability and H-bonding of water, provides an excellent description of the solvation and dynamics of hydroxide and hydronium in bulk water. Here we turn our attention to slabs, cylinders, and droplets. In extended simulations with 1000 molecules, we find that hydroxide consistently prefers the surface, hydronium consistently avoids the surface, and the two together form an electrical double layer until neutralization occurs. The behavior of hydroxide can largely be accounted for by the observation that hydroxide moving to the surface loses fewer hydrogen bonds than are gained by the water molecule that it displaces from the surface. At the same time, since the orientation of the hydroxide increases the ratio of dangling hydrogens to dangling lone pairs, the proton activity of the exposed surface may be increased, rather than decreased. Hydroxide also moves more rapidly in the surface than in the bulk, likely because the proton donating propensity of neighboring water molecules is focused on the one hydrogen that is not dangling from the surface.

Introduction

Water molecules occasionally autoionize, forming hydronium and hydroxide ions with surface propensities that are important in a wide range of contexts. However, despite intensive experimental and theoretical efforts, it has proved vexing to establish the surface population of ions, with some results suggesting excess hydroxide and others excess hydronium.[1−5]

It has long since been known that air bubbles and oil droplets in water respond to electric fields as though negatively charged, and experimentalists have gone to great lengths to exclude artifacts due to impurities.[6−10] The simplest explanation of this electrophoretic behavior is that hydroxide accumulates at the interface while hydronium is located on the other side of the slip plane. This interpretation is supported by the time constant for the decline of the surface tension from the value for a newly prepared surface to the equilibrium value.[11] It is also consistent with the variation of the surface tension with pH and salt concentration.[12−14]

Surface-sensitive spectroscopy has also been used to study the surface character of water, including second harmonic generation (SHG),[15−17] vibrational sum-frequency spectroscopy (VSFS),[18] and photoemission spectroscopy (PES).[19] However, interpretation of data from those techniques is difficult. One problem is uncertainty about the range of depths being probed.[1] Another concern is that detection of the spectroscopic signature of hydroxide requires high concentrations achievable only by the addition of counterions[1,20] at concentrations that are likely to be significantly perturbative.[11,21] An analogous concern arises from the high surface concentration of pH indicator required for spectroscopic measurement of surface pH.[22]

Theoretical results have proven highly sensitive to the approach taken. First-principles calculations depend on the choice of density functional: different functionals give very different pictures of ion solvation and ion dynamics.[23] First-principles calculations also impose severe practical limitations on the numbers of molecules, concentrations of species, and durations of simulations. Previous DFT studies provide mixed results for the air–water interface. Some of them predict a preference of hydronium for the surface,[24,25] some show a slight hydroxide enhancement at the surface,[26] and some conclude that there is no significant surface preference for either ion.[27]

Molecular mechanics and MS-EVB results also depend on the model used[28] and the parametrization[29] of hydroxide vs hydronium vs water. Several studies using MS-EVB have hydronium favoring the surface[30−33] and hydroxide avoiding it.[34] However, in a recent study, Wick showed that hydronium’s attraction to the surface in previous MS-EVB studies was due to neglect of polarizability.[28] Moreover, Wick and Dang showed that their polarizable MS-EVB potential predicts little free-energy change in moving hydroxide to the surface.[35]

Recently, water has been modeled with the LEWIS force field in which valence electron pairs are efficiently described as explicit semiclassical particles that interact pairwise with each other and with semiclassical kernels that comprise nuclei and core electrons. The resulting water is amphiprotic and intrinsically polarizable.[36] No reparametrization is involved in describing different protonation states or intermediates between them. A study of the self-ions shows hydroxide solvation that largely mirrors that of hydronium (with the Eigen form most common and the Zundel form as the transition state for proton transfer) and ion diffusion rates relative to water that are consistent with experiment.[37]

Here we use LEWIS to simulate systems with an ion among 1000 molecules in bulk (B), slab (S), cylinder (C), and droplet (D) configurations (i.e., volume concentration ∼0.05 M or surface concentration ∼1 ion/10 nm2). To avoid bias in positioning the ions in systems with broken symmetry, we condense from the gas state at 1000 K to a slab, cylinder, or droplet at 300 K. For each box, we first equilibrate at 1000 K for 500 ps to generate a fully distributed gas state. A subsequent shift to 300 K resulted in a slab, cylinder, or droplet. Ten systems were studied as follows:

protonated to result in 1 hydronium and 999 water molecules.

deprotonated to result in 1 hydroxide and 999 water molecules.

1000 neutral water molecules.

protonated to result in 1 hydronium and 999 water molecules.

deprotonated to result in 1 hydroxide and 999 water molecules.

1 hydronium, 1 hydroxide and 998 water molecules.

protonated to result in 1 hydronium and 999 water molecules.

deprotonated to result in 1 hydroxide and 999 water molecules.

protonated to result in 1 hydronium and 999 water molecules.

deprotonated to result in 1 hydroxide and 999 water molecules.

System S0 was designed to investigate the ion-free process of condensation from a 1000 K gas to a 300 K slab and to characterize H-bonding within the slab. Systems S+ and S– were used to study the distribution of hydronium and hydroxide ions, respectively. We used system S± to compare the behavior of coexisting ions with the behavior of the solitary ions in systems S+ and S–. Systems B+ and B– were built to compare the H-bonding of ions in the bulk with that in the slab. We also constructed systems C+, C–, D+, and D– to test the influence of surface curvature.

Results and Discussion

Figure a shows a typical condensation in system S0. Before condensation, water molecules are distributed throughout the box. When the temperature is decreased, water molecules begin to form clusters and concentrate in one region of the box (translated to the middle after the simulation). After 60 ps a stable slab is formed. Figure b shows the evolution of the number of H-bonds and the distribution of oxygen atoms. The rmsd is for the z-coordinates of oxygen atoms and decreases from the value expected for a uniform gas to that expected for a uniform slab. Meanwhile, the number of H-bonds increases from a small number due to some dimers in the gas to 1.5 per water molecule in the condensed phase.

Figure 1

Condensation in system S0. (a) Configurations starting from the homogeneous gas state previously equilibrated at 1000 K. Oxygen atoms are red, protons white, and electron pairs green. (b) Number of H-bonds (black) and root-mean-square deviation along the long axis of the box (red) for the condensation shown in panel a. Circles at t = 0 in panel b indicate starting values. Final values are as expected for a homogeneous slab at the density of bulk water.

Ion Distribution

To obtain initial statistics for the surface propensities of the ions, nine parallel simulations were performed for systems S+ and S–. Figure a–b shows the z-coordinates of the oxygen in each ion relative to the average for all the molecules (i.e., the middle of the slab). At t = 0 ps, before condensation starts, ions could appear at any position in the evenly distributed gas state. After cooling, ions move with water molecules as the slab begins to form. Figure a shows that, after 150 ps, hydronium ions are distributed at many depths, but avoid the surface (located at z – ⟨z⟩ ∼ ±15.5 Å). In contrast, Figure b shows that, after 150 ps, all hydroxide ions are located at the slab surface, independent of their earlier positions. These preferences did not change for the coexisting anion and cation in system S± (Figure c). Until they neutralize (not until ∼300 ps in this particular simulation), the ions form a double layer.

Figure 2

Ion distributions relative to the middle of the slab. Trajectories for (a) hydronium in system S+, (b) hydroxide in system S–, and (c) both ions in system S±. The 9 condensations for systems S+ and S– were each started from a different configuration of the gas phase that had been fully equilibrated at 1000 K. The single trajectory for S± was obtained by removing a proton from an equilibrated S+ slab to create a buried hydroxide that is seen to rapidly move to the surface and stay there. (d) Statistics for water (black), hydrogen bonds (purple), hydroxide (red), and hydronium (blue). Water and hydrogen bond statistics were derived from 5000 snapshots obtained at 10 fs intervals from the final 50 ps of the simulation of system S0. The vertical dashed line identifies the Gibbs dividing surface. Hydroxide statistics (red) were derived from 4500 snapshots obtained at 100 fs intervals from the final 50 ps of 9 parallel simulations of system S–. Hydronium statistics (blue) were obtained from 29,500 snapshots obtained at 100 fs intervals from the final 50 ps of 57 parallel simulations of system S+. The bin width is 0.1 Å.

Although the hydronium ions in our initial nine trajectories seem to avoid the surface, as shown in Figure a, this could be a coincidence. To obtain better statistics, we performed 48 additional simulations for system S+ and collected data from the last 50 ps of each trajectory. These results appear in Figure d, along with those for water and hydroxide. Here, the black curve shows the uniform distribution of water molecules, through the interior of the slab, with the decline beyond 14.5 Å defining the surface. As expected, the Gibbs dividing surface (vertical dashed line) is located at 15.6 Å. The purple curve shows the variation of the ratio of the average number of H-bonds in each bin to the average number of oxygens in that bin. This measure also shows that the interior of the slab is uniform, although H-bonding is somewhat disrupted in the water layer just inside the surface. In sharp contrast, the red curve shows that hydroxide strongly prefers the surface and the blue curve shows that hydronium strongly avoids the surface. The variations of the blue curve at greater depths are more difficult to interpret. While our slab contains only ∼5 water layers from each surface, the water and H-bond distributions cited above suggest that the surface does not perturb the structure of the water beyond the first layer inside the surface. However, perturbations by hydronium may depend on distance from the surface. Another possibility is that 57 parallel trajectories are still not enough to generate a flat distribution for hydronium.

Influence of an Electric Field

The surface adsorption energy for hydroxide has been estimated to be 16–25 kBT, from the pH dependence of the ξ potential of oil droplets in water[38] and from an analysis of the relaxation of the surface tension of the air–water surface.[11] These values are impressive since they are larger than the 40.65 kJ/mol heat of vaporization of water,[39] and presumably represent significant effects of charge.

To test the surface adsorption energies of both ions using the LEWIS model, we applied an electric field to systems S+ and S–. We chose 0.5 V/nm because, at 300 K, 20 kBT amounts to 4.98 × 104 (V·C/mol) or 0.5 V per electronic charge (9.65 × 104 C/mol). Our chosen electric field provides that energy difference in 1 nm.

We performed simulations with 0.5 V/nm electric fields oriented in the positive and negative z directions, starting from the last frame of well-equilibrated slabs in S+ and S– systems. Statistics for the rmsd and the number of H-bonds during these simulations (not shown) indicated that the electric fields do not have a significant effect on the density of the slab or the structure of the H-bond network.

Figure a shows the hydronium position under the influence of electric fields in opposite directions. We see that 0.5 V/nm suffices to drive the ion across the slab, but not to drive it into the surface region on either side. On the other hand, Figure b shows that 0.5 V/nm can shift the hydroxide slightly within the slab surface, but not drive it out of the surface. These results indicate that the barrier to hydronium moving into the slab or hydroxide moving out of the slab is steeper than 20 kBT/nm. We also tried doubling the E-fields to 1.0 V/nm. But such strong fields caused the slabs to break down. In those simulations, hydronium remains solvated in water clusters while hydroxide escapes and travels across the box boundaries.

Figure 3

Ion trajectories under the influence of electric fields perpendicular to the slab in each direction: (a) hydronium in system S+ and (b) hydroxide in system S–. The magnitude of the electric field in each case is 0.5 V/nm.

Driving Force

To gain insight into the molecular basis for the surface propensity of both ions, we begin by comparing the H-bonding of all three species in the bulk and at the surface. Based on the water distribution in Figure d, it is reasonable to define the outmost layer |z – ⟨z⟩| > 14.5 Å as the slab surface in systems S+ and S–. For bulk statistics we use trajectories for systems B+ and B–. The results are shown in Figure .

Figure 4

The average number of H-bonds accepted (gray) and donated (black, imperceptible for hydroxide) by the indicated species, in the slab surface (S) vs in the bulk (B). The slab statistics were obtained from the same trajectories as used in Figure d. We consider |z – ⟨z⟩| > 14.5 Å as the surface region in the slab. The bulk statistics are obtained from systems B+ and B–.

In principle, possibilities for donating:accepting H-bonds in bulk water are 2:2 for H2O, 3:1 for H3O+, and 1:3 for OH–. The balance in H2O is expected to lend uniformity and strength to the H-bond network, and the imbalance H3O+ and OH– is expected to be disruptive. In practice, we find that all three species in the bulk participate in only three H-bonds, on average, at 300 K. In the case of H2O, donating and accepting remain well balanced at ∼1.5:1.5. However, the skew increases to ∼3:0 for H3O+ and ∼0:3 for OH–. These considerations alone suggest that both H3O+ and OH– might be better situated at the surface, especially if they are able to each form H-bonds with three neighboring waters molecules and if the water molecules that they displace from the surface are able to form more H-bonds in the bulk.

The surface statistics in Figure show that both water and hydroxide gain some H-bonds in moving from the surface to the bulk. However, the extent differs dramatically: on average, while water gains 0.8690 H-bond, hydroxide gains only 0.2805. These numbers indicate that, when a hydroxide displaces a water molecule from the surface, there is a net gain of 0.5885 hydrogen bond. This is qualitatively consistent with the tendency of hydroxide to prefer the surface.

For a more quantitative assessment, we note that the heat of vaporization of water is 40.65 kJ/mol.[39] Dividing this number by an average of 1.5 H-bonds/water molecule (see Figure ) we estimate that each H-bond costs about 27.1 kJ/mol, such that, when hydroxide displaces water from the surface, the energy decreases by 16.0 kJ/mol. This value is smaller than the about 20 kBT adsorption energy estimated by experiments.[11,38] The underestimate may be due to limiting our analysis to H-bonds inside the first solvation shell. Nevertheless, at 300 K, where kBT = 2.49 kJ/mol, our estimate of 16.0 kJ/mol gives a dramatic surface/bulk Boltzmann factor of 590 for hydroxide ion.

Unfortunately, we are unable to carry out a similar analysis for hydronium ions because they are never found at the surface in our simulations. However, their absence suggests that, whereas surface hydroxide is able to accept almost three H-bonds, surface hydronium is unable to accept an adequate number of H-bonds, contrary to suggestions elsewhere.[17,33] This may be due to the difference between the diffuse, ringlike distribution of the lone pair electrons in hydroxide vs the relatively rigid, tetrahedrally directed distribution of the protons in hydronium.

Hydroxide Orientation and Surface pH

The reactivity of surface species should be related to their orientations. As shown in Figure , although water loses H-bonds at the surface relative to the bulk, it loses less than half of them and it still donates and accepts H-bonds in approximately equal numbers, such that the numbers of “dangling” lone pairs are similar to the numbers of “dangling” protons. This balance between donated and accepted H-bonds is consistent with a surface H-bond network that lends strength to the surface tension of water and contradicts DFT simulations that indicate large differences between numbers of donated and accepted H-bonds in the vicinity of the surface.[40]

In the case of hydroxide, very few H-bonds are lost at the surface relative to the bulk. The ion continues to accept an average of almost three H-bonds and must therefore be oriented with all lone pairs buried. Figure shows statistics for θ, defined as the angle between the O–H vector and the slab normal, as pictured in the figure inset. It is clear that the hydrogen atom of the hydroxide always points to the air side. The distribution is well fit bywhich assumes that the energy increases quadratically with θ and takes the sin(θ) dependence of the cone size into account.

Figure 5

Distribution of hydroxide ion orientations. The angle θ is defined in the inset. Bin size is 1°. Statistics were obtained from the same 9 trajectories of system S– as used in Figure d. The solid line shows the fit to eq .

The irony is that, while hydroxide’s preference for the surface gives the surface a negative charge, the ion’s orientation in the surface renders the outermost surface somewhat enriched in dangling hydrogens and somewhat depleted in proton acceptors (i.e., dangling lone pairs). Thus, since pH properly refers to proton activity, water may present a mildly acidic surface even though it harbors excess hydroxide.

Hydroxide Surface Mobility

Much attention has been paid to the diffusion of water ions in the bulk liquid. Whereas the diffusion of water requires movement of the whole molecule, its ions can diffuse by proton hopping. The diffusion of hydronium is the fastest, as it requires only that the excess proton move from one water to a neighboring water. For hydroxide, proton hopping is more difficult as it involves breaking up a neighboring water molecule. Figure shows that hydroxide diffuses more rapidly on the surface than it does in the bulk. Whereas D = (1/6)d⟨r2⟩/dt = 0.15 Å2 ps–1, in the bulk, on the surface D = (1/4)d⟨r2⟩/dt = 0.28 Å2 ps–1. This is consistent with the fact that water molecules in the surface are each donating just one hydrogen bond. With one proton dangling from the surface, the propensity to donate is focused entirely on the other proton, making it more readily available to transfer to a hydroxide.

Figure 6

Diffusion of hydroxide. The mean square distance vs time is shown for hydroxide in the bulk (solid line) and hydroxide at the surface (dashed line). The latter was derived from the last 125 ps of the 7 trajectories in Figure b in which hydroxide arrives at the surface relatively quickly.

Curvature Dependence

As has generally been the case elsewhere, the above results concern the overall flat air–water interface of a slab, and it is worth considering whether ion behavior might be different at a curved interface. In order to evaluate the influence of curvature we varied box dimensions to condense water into cylinders and droplets that contain one ion (systems C+, C–, D+, and D−).[41] As for the slabs, systems were equilibrated at 1000 K for 500 ps and then condensed and equilibrated at 300 K for another 500 ps. The linear and parabolic water distributions shown in black in Figures a and 7b are characteristic of cylinder and droplet shapes, respectively. In both the cylinder and the droplet, as in the slab, hydroxide prefers the surface and hydronium avoids it.

Figure 7

Distribution of ions relative to curved surfaces. (a) Ion distances from the central axis of a cylinder. (b) Ion distances from the center of a droplet. Color code is as in Figure d. Statistics were obtained from snapshots at 100 fs intervals during the last 100 ps of five parallel 500 ps simulations of systems C+, C–, D+, and D–. The bin width is 0.1 Å.

Methods

In this work, boxes with periodic boundary conditions were constructed for studying 1000 molecules in bulk (B), slab (S), cylinder (C), and droplet (D) configurations. The dimensions wereNote that the density of 1000 water molecules in a box of type B is 997.05 g/cm, i.e, the experimental density at 300 K and 1 atm.[39]

We use the reactive and polarizable LEWIS model of water, comprising independently mobile valence electron pairs, protons and oxygen kernels, with a fictitious electron mass of 1 amu and interactions truncated with compensation at 9 Å.[36] The pairwise potentials describing the interactions between these particles were trained on the structures and relative energies of neutral, protonated, and deprotonated water monomers and dimers.[42] The resulting force field provides an excellent description of Grotthuss transport in water chains,[42] the polarization of water in the neat liquid,[36] and the correct trend for diffusion of water, proton holes, and excess protons in bulk water.[37]

The NVT molecular dynamics simulations were performed in the GROMACS package (version 4.6.3)[43] using the velocity-rescaling thermostat[44] with a time constant of 0.01 ps. The time step was 0.2 fs for all simulations. Equilibration for 500 ps was at 300 K for the bulk systems and at 1000 K for systems to be condensed into anisotropic systems. Cooling of the latter to 300 K generally resulted in a stable slab after ∼60 ps, and a stable cylinder or droplet after ∼350 ps. The delay in the latter systems is due to the larger sizes of the boxes. For S0, snapshots were stored every 10 fs over 500 ps. For systems S+, S–, S±, B+, and B–, snapshots were stored every 100 fs over 200 ps trajectories. For systems C+, C–, D+, and D–, snapshots were stored every 100 fs over 500 ps trajectories.

Counts of numbers of H-bonds depend on the water model and the somewhat arbitrary definition of an H-bond.[45] Using an H-bond defined by an attraction of 9.4 kJ/mol or more between two water molecules, Jorgenson et al. obtained an average of 3.59 H-bonds for a TIP4P water molecule at 298 K.[46] More recently, defining an H-bond by (1) a H···O distance no longer than 0.25 nm (i.e., the first minimum on the O–H radial distribution function) and (2) an O–H···O angle exceeding 150°, Zielkiewicz obtained an average of 2.357 H-bonds for a TIP4P water, 2.457 for a SPC/E water, and 2.082 for a TIP5P water, at 298 K.[47] In this work, we use the GROMACS[48] requirements of (1) an O···O distance smaller than 0.35 nm and (2) a H–O···O angle not larger than 30°. These criteria gave an average of 3.0 H-bonds for a LEWIS water molecule, as shown in Figure .