Characterizing the bifurcating configuration of hydrogen bonding network in interfacial liquid water and its adhesion on solid surfaces.

The interfacial structures of liquid water molecules adjacent to a solid surface contribute significantly to the interfacial properties of aqueous solutions, and are of prime importance in a wide spectrum of applications. In this work, we use molecular dynamics (MD) simulations to explore the interfacial structures, mainly in term of hydrogen bonding network, of a liquid water film interacting intimately with solid surfaces, which are composed of [100] face centered cubic (FCC) lattices. We disclose the formation of a bifurcating configuration of hydrogen bonds in interfacial liquid water and ascribe its occurrence to the collective effects of water density depletion, hydrogen bonds and local polarization. Such bifurcating configuration of interfacial water molecules consists of repetitive layer by layer water sheets with intra-layer hydrogen bonding network being formed in each layer, and inter-layer defects, i.e., hydrogen bonds formed between two neighboring layers of interfacial water. A lower bound of 2.475 for the average number of hydrogen bonds per interfacial water molecule is expected. Our MD study on the interfacial configuration of water on solid surfaces reveals a quadratic dependence of adhesion on the solid-liquid affinity, bridging the gap between the macroscopic interfacial property W adh and the microscopic parameter ε SL of the depth of the Lennard-Jones solid-liquid potential. This journal is © The Royal Society of Chemistry.

Introduction

The wettability of water or aqueous solutions on solid surfaces is crucially dependent on the interfacial properties of water interacting intimately with the surface materials. In the past decades, rapid developments in nanomaterials and microfabrication create an increasing demand for insightful understanding of both the static and the dynamic behaviors of liquid water near solid surfaces, due to its critical importance in a wide spectrum of applications, such as surface coating,[1,2] micro/nano-fluidic devices,[3-5] phase change heat transfer,[6,7] and various biological occasions.[8] For liquid water residing on a solid surface, the adhesive force applied by the solid surface, as well as the cohesive character of water, induces the reconfiguration of water molecules adjacent to the solid surface and in a large part determines the interfacial properties.[9] For example, the parallel dipole alignment of water molecules near solid surfaces has been demonstrated to minimize the loss of hydrogen bonds and in the meantime contribute substantially to the interfacial tension.[10]

Depending on the strength of adhesion, solid surfaces can be categorized into two main groups: hydrophobic and hydrophilic.[11,12] A facile way to reflect the surface wettability is via its contact angle derived from the Young's equation: , where γ is the interfacial tension and the subscripts of S, L and V represent solid, liquid and vapor phases, respectively. However, this equation does not provide any insightful information on the microscopic details near the solid–liquid interface, especially the structural change of liquid water in response to surface wettability.[13,14] Increasing interests in micro/nanofluidics have inspired studies on solid–liquid interfaces whose characteristic length drops down to nanometers. First and foremost, a depletion region of liquid water is formed near a solid surface and such density deficit becomes more prominent for hydrophobic surfaces.[15,16] Above the depletion region, wave-like spatial undulations of liquid density have been reported,[17-19] implying a new mechanism of energy dissipation in the form of layer-by-layer friction during the process of dynamic wetting.[20-22] It was also reported that such spatial undulations can break the ideal tetrahedron geometry of hydrogen bonds,[23] which is usually formed in bulk water, and may give rise to the polarity of interfacial water.[24,25] In this respect, the identification of truly interfacial molecules (ITIM) analyses revealed that the effect of structural change on interfacial properties may be limited to the close vicinity of the interface, i.e., the first molecular layer.[26,27] Furthermore, when electrostatic interactions are taken into account, the balance or symmetry of accepting and donating hydrogen bonds is broken in the interfacial layer, leading to surface charging.[28]

In this study, we use molecular dynamics (MD) simulations to explore the interfacial structures and in particular the hydrogen bonding network of a liquid water film interacting with solid surfaces. We show that the bifurcating configuration of interfacial water, which is jointly determined by water density depletion, hydrogen bonds and local polarization, contribute significantly to the interfacial properties of liquid water. Importantly, this MD study reveals a quadratic dependence of work of adhesion Wadh, a macroscopic and experimentally accessible value, on the microscopic parameter εSL of the depth of the Lennard-Jones solid–liquid potential, indicative of solid–liquid affinity. Our findings can advance the fundamental understanding on the properties of solid–water interface.

Materials and methods

In our MD simulations, a water film of 10 nm wide and 5 nm thick was placed on a solid surface. Water molecules were modeled as a rigid planar four-site system, i.e., the TIP4P water model,[29,30] which consists of the 12-6 Lennard-Jones interaction among oxygen atoms with the energy parameter εOO and the distance parameter σOO, as well as partial charges on oxygen and hydrogen atoms. We also tested the influences of water models by comparing the results based on the TIP4P model to those based on the SPC/E model (ESI†), and the results of both models are generally identical. The solid surface consisted of two layers of [100] face centered cubic (FCC) lattices and a refined Lennard-Jones model (εSS, σSS), which is able to accurately reproduce gold–water interface, was used to model the interactions among solid atoms.[31] During the simulations, solid atoms were tethered to their initial positions via a harmonic potential, such that their mean displacement was small enough to mimic a rigid wall.[32] The first layer of solid atoms in contact with the water film were slightly randomized in the surface normal direction with an amplitude of 0.3σSS to avoid the unrealistic large slip length observed in previous MD simulations.[21,33] At first, the solid–water interactions were modeled via the 12-6 Lennard-Jones potential (εSL, σSL) with geometrical combination rule, i.e., εSL = 3.8 kJ mol−1 and σSL = 0.298 nm. In order to study the structural response of water molecules to surface wettability, we tuned εSL from 3.8 kJ mol−1 to 1.0 kJ mol−1 in this study. It has been shown that decreasing εSL can effectively enhance surface hydrophobicity and this method has been widely used in previous studies.[34,35]

All of our MD simulations were carried out on Gromacs 5.1.2.[36] Provided that the realistic density of liquid water is 996.5 kg m−3 at 300 K, there were 15 300 water molecules included in the modeling domain. The time step was set to be 1 femtosecond (fs). Short range interactions were cut-off at 1.2 nm with the Verlet scheme. The particle mesh Ewald method with cubic interpolation and grid spacing of 0.12 nm was adopted to handle the long range Coulombic interaction.[37] An NVT ensemble with Nose–Hoover thermostat was used to maintain the system temperature at 300 K. Firstly, the system was equilibrated for 10 ns and then the simulation was continued for another 10 ns for data collection. In our simulations, the criteria for determining a hydrogen bond are based on the following geometrical conditions:[38]

(1) The distance between two oxygen atoms is not greater than RcOO = 0.34 nm;

(2) The distance between the acceptor oxygen and the donor hydrogen is not greater than RcOH = 0.2425 nm;

(3) The H–O⋯O angle is less than αcH−O⋯O = 30°.

Results and discussion

Fig. 1(a) presents the density profiles of oxygen atoms for interfacial water molecules in response to varied wettability. The density of oxygen atoms ρ is normalized by its value ρ0 in bulk water. For a liquid–vapor interface, water density shows a monotonic increase from zero in the vapor phase to its bulk value in the liquid phase. Thus-formed density profile can be accurately predicted by a hyperbolic tangent function indicative of a diffusion-governed mass transport.[39,40] However, the density distribution at the solid–liquid interface exhibits a non-monotonic configuration, in spite of different εSL. The local density rapidly increases to a peak value, then oscillates around the bulk value ρ0 with decaying amplitude and eventually settles to ρ/ρ0 = 1. The strong density undulations in Fig. 1(a) provide a clear evidence of near-surface water layering, i.e., water molecules are mainly trapped within each spike in the density profile. Such layering effect extends about 1 nm into the water film. Within this range, the resultant layer-by-layer structure implies a mass transport mechanism jointly governed by phase transition,[20,41] convective mass transfer[42] and self-diffusion.[43] In addition, lowering εSL mitigates the water layering and may lead to a smooth transition to a liquid–vapor like interface when the solid surface becomes more hydrophobic.

Fig. 1

(a) Density of oxygen atoms ρ normalized by the bulk density ρ0 as a function of the distance z from the solid substrate. The dashed line is the density profile of water at the liquid–vapor interface. (b) Distribution of number of hydrogen bonds per water molecule, nHB(z), as a function of the distance z from the solid substrate.

We then estimated nHB, the number of hydrogen bonds shared by each water molecule, as a function of the distance z from the solid surface with different εSL. In bulk water, water molecules at equilibrium are apt to retain a local ice-like tetrahedral structure and each water molecule is expected to form 4 hydrogen bonds with its adjoining neighbors. In reality, nHB is smaller than 4 and varies from 3.2 to 3.5, since hydrogen bonds in liquid water are constantly breaking and reforming.[10,19] In our MD simulations, nHB is found to be 3.3 on average in the bulk water. Fig. 1(b) demonstrates that liquid water experiences a loss of hydrogen bonds when it moves closer to the solid surface and increasing εSL can effectively reduce the loss of hydrogen bonds. For the standard case of εSL = 3.8 kJ mol−1 (strong hydrophilicity), the number of hydrogen bonds is almost unchanged comparing to that in the bulk water.

Hydrogen bonding in water is in fact a coulombic attraction between partially charged oxygen and hydrogen atoms. Therefore, the mathematical expectation of nHB for a water molecule can be summarized as:where p0 is the proportionality constant, ρO1 and ρO2 are the normalized density of oxygen atoms within RcOH from two hydrogen atoms, and ρH is the normalized density of hydrogen atoms within RcOH from the oxygen atom. A rough estimation for p0 is 0.825, since a water molecule is expected to form 3.3 hydrogen bonds on average with its neighbors in bulk water. Essentially, the water layering structure as shown in Fig. 1(a) implies that a water molecule is most likely to donate its hydrogen atoms to accepting oxygen atoms from the same layer. In this scenario, water molecules within the same spike are preferably connected by intra-layer hydrogen bonds, forming a planar sheet.[44]

To testify this hypothesis, we examined the height difference Δz between the donating and accepting oxygen atoms for each pair of water molecules connected by hydrogen bonds near the solid surface. Hereafter, hydrogen bonds are considered to be near the solid surface if those connected water molecules are located within the first two layers indicated by Fig. 1(a). Fig. 2(a) presents the probability density distribution of Δz ranging from −RcOO to RcOO. In bulk water, the probability density profile of Δz is uniform, suggesting that a water molecule is equally probable to form hydrogen bonds with peers from all directions. In contrast, the probability density of Δz near the solid surface exhibits a pulse-like profile peaking at Δz = 0, indicating that water molecules in each spike of Fig. 1(a) are preferably aligned in parallel to the solid surface and form a planar sheet. In addition, increasing surface wettability intensifies the layering effect and subsequently stabilizes the planar configuration of each water layer, as shown in Fig. 2(a). To provide a detailed view of the as-formed planar water sheet, we analyzed the triple angle β, which is defined as the angle formed between two neighboring oxygen atoms connected by intra-layer hydrogen bonds (illustrated in the inlet of Fig. 2(b)). Fig. 2(b) presents the probability density distribution of β and highlights two most probable values at β = 90° or 180°. Therefore, the planar network of interfacial water molecules can be represented as a water sheet comprised of square cells. Although the layering structure of interfacial water has been widely observed, the planar geometry of each water layer (sheet) is not unambiguously a square grid system and may vary depending on the lattice structures of the solid surface. For example, a hexagonal configuration rather than a square grid is expected for interfacial water molecules close to a graphene surface.[45]

Fig. 2

(a) Probability density profiles of the distance Δz (Δz ranging from −RcOO to RcOO) between the accepting and the donating oxygen atoms along the outer normal direction of solid surfaces. (b) Probability density profiles for the triple angle β, which is defined as the angle between the two oxygen atoms connected by hydrogen bonds (stood by the red dashed bars in the inlet).

Despite the stabilization effect of strong confinement from the solid surface atoms, as-formed water sheets with square cells can still be distorted by viscous interactions among water molecules, mostly in the form of hydrogen bonds. A secondary configuration of interfacial water implied by Fig. 1(a) is hydrogen bonds bridging two neighboring water sheets. Since the distance between neighboring density peaks is around 0.28 nm (smaller than RcOO = 0.34 nm) in Fig. 1(a), the donating hydrogen atoms can be attracted by the accepting oxygen atoms from the adjacent layers. In this case, inter-layer hydrogen bonds may probably form, disturbing the otherwise planar geometry of each water sheet. In this respect, we name such secondary configuration of hydrogen bond as inter-layer defect. Those inter-layer defects are manifest in Fig. 2(a) as the two minor peaks located at Δz = ±0.28 nm and become more prominent with increasing εSL. Although such inter-layer defects are seemingly consistent with those perpendicular hydrogen bonds proposed by Velasco-Velez et al.,[44] the asymmetric profiles in Fig. 2(a) reveal a directional preference of as-formed inter-layer defects. Since Δz is more likely to be positive, the inter-layer defects tend to connect donors within the current layer to acceptors from the subsequent layer.

We show that the directional preference of inter-layer defects can be explained by the local polarization of interfacial water. A force analysis on a water molecule interacting with the solid surface is illustrated in Fig. 3(a). FH and FO denote the forces that the gold substrate applies to the oxygen and hydrogen atoms of a water molecule, respectively. In most cases, FH is negligible comparing to FO and is taken to be zero in this study. Therefore, the interaction between a water molecule and the solid substrate can be simplified as a force FO acting at the center of mass (COM) of the water molecule together with a resultant torque T⃑O = F⃑O × L⃑OC, where L⃑OC is the vector connecting the oxygen atom and the COM. Here, FO can be evaluated as:where ρS is the number density of solid (gold) atoms and Rc is the cut-off radius. From eqn (2), FO can be positive (repulsive), or negative (attractive), depending on the distance z of the oxygen atom from the solid surface. When , FO is positive and water molecules within this region is hence repulsed away, leading to a depletion region with a thickness of z*. In this scenario, the water film is seemingly suspending z* from the gold surface. As FO is attractive for z > z*, T⃑O will force the water molecule in Fig. 3(a) to spin in the counterclockwise direction and subsequently induces the local polarization of interfacial water.

Fig. 3

(a) Schematic illustration of the orientation of a water molecule. αHH stands for the angle between the H–H bond vector and the surface normal. αOH indicates the angle made by the O–H vector with the surface normal. (b) Probability density profiles of H–H bond vector. (c) Probability density profiles of OH vector. (d) Illustration of bifurcating hydrogen bonding network near solid–water interface.

We examined the orientation of a water molecule close to the solid surface in terms of two vectors, i.e., the H–H bond vector, which points to the opposite direction of its dipole moment, and the O–H vector, which connects the oxygen atom to one of the hydrogen atoms that are bonded to it. In this way, the orientation of a water molecule can be determined by two angles, i.e., αHH between the H–H bond vector and z axis, and αOH between the O–H vector and the z axis. Fig. 3(b) and (c) present probability density profiles of αHH and αOH for different surface wettability. Without the strong confining effect of the solid surface, the probability density distributions for both αHH and αOH ought to be homogenous, since a water molecule shows no orientational bias in the bulk. However, the highly distorted profiles in Fig. 3(b) and (c) provide unambiguous evidences for the orientational change of interfacial water molecules.[44] The central peaks imply that interfacial water molecules are mostly aligned in parallel to the solid surface and hence form a planar sheet. In the meantime, a closer look at Fig. 3(b) and (c) reveals a preferential orientation of αHH, αOH ≤ 90°, which agrees with our force analysis on the local polarization. Moreover, such polarization favors the upright orientation of water molecules, and is responsible for the directional preference of inter-layer defects. We also notice that there are irregular wiggles that fall within the left sides of Fig. 3(b) and (c). These wiggles produce a secondary peak of the probability density profile, which corresponds to the inter-layer defects. Considering the H–O–H bond angle αHOH within a TIP4P water molecule is taken as 104.52°, those secondary peaks occur right at and αOH = αHOH − 90°. Therefore, a typical configuration of inter-layer defects and intra-layer network is illustrated in Fig. 3(d).

So far, we perceive a bifurcating configuration of interfacial water, i.e., lateral water sheets and perpendicular inter-layer defects, as shown in Fig. 3(d), due to the collective effects of water density depletion, hydrogen bonds and local polarization. According to eqn (2), the first layer of water molecules is formed in such a fashion that the strong repulsive force leads to a depletion region and water molecules above it are attracted toward the upper edge of the depletion zone. The aggregation of water molecules in this layer produces a cloud of oxygen atoms, which drives the originally randomly oriented O–H bonds to be aligned laterally so as to form intra-layer bonds. Eventually, the aggregated water molecules in each spike are fashioned into a planar water sheet. For the [100] FCC lattice structures of the solid surface used in this work, the water sheet is ideally comprised of square-like lattices and one water molecule is expected to form 4 (2 donating and 2 accepting) hydrogen bonds with its in-plane neighbors. However, the local polarization may break one of the donating hydrogen bonds and expel one hydrogen atom to the upper space. This dispersed hydrogen atom may or may not form a hydrogen bond with subsequent water molecules. Given that in the interfacial layer, a lower bound for the average number of hydrogen bonds per water molecule can be estimated as:

Eqn (3) implies that the contribution of hydrogen bonds to the solid–liquid interfacial tension γSL is relatively trivial comparing with that to the liquid–vapor interface. Following the first layer formation, a second layer of water molecules comes into being once the inter-layer hydrogen bonds are formed between the dispersed hydrogen atoms and the subsequent water molecules. Similar to the first layer, the second layer is connected by intra-layer hydrogen bonds and inextricably expel hydrogen atoms into the upper space as well, giving rise to the third layer formation and vice versa. Therefore, a repetitive pattern of water sheets straddled by inter-layer defects is formed. Due to the fact that the polarization effect T⃑O decays rapidly as T⃑O ∼ 1/z4, the amplitudes of the resultant density peaks in density profiles diminish successively and gradually settle at ρ = ρ0.

We hereby propose a damped harmonic oscillator function to describe the non-monotonic density profiles of interfacial water and quantify the bifurcating configuration of hydrogen bonds as:where A is the amplitude, B the damping coefficient, λ the inter-layer spacing, and δ the depletion length. As shown in Fig. 4(a), the proposed function yields a satisfactory agreement with the density fluctuations near the solid surface. Among unknown parameters, λ is found to be 0.28 nm for TIP4P water by analyzing the optimal distance between oxygen atoms that are connected by hydrogen bonds in bulk water (ESI†). Note that λ = 0.28 nm is also in line with the reported location of secondary peaks in Fig. 2(a). The depletion length δ is defined as the upper edge of the depletion zone and corresponds to the location of first density peak. If δ is solved by taking FO = 0 in eqn (2), a constant value of is expected for all cases. However, it is found that the depletion length decreases with increasing εSL, as shown in Fig. 4(a). This is because the cohesive force FLL of water must be taken into consideration in determining δ. Since FLL arises from the interactions among water molecules, it can be taken as constant regardless of εSL. Therefore δ can be solved by equating the right hand side of eqn (2) to FLL:

Fig. 4

(a) Density profiles fitted to . (b) Work of adhesion Wadh and depletion length δ with respect to εSL. The quadratic increase of Wadh with εSL is attributed to εSL ∼ δ−2.

Apparently a smaller δ is expected for increasing εSL. Note that eqn (5) can not be solved explicitly. Therefore, by pinpointing the locations of first density peaks associated with different values of εSL, a simple scaling law relating δ to εSL as shown in Fig. 4(b) is unveiled:

As-formed bifurcating configuration also has a profound influence on the interfacial properties of water. We consider the effect of such interfacial molecular structures by evaluating work of adhesion as:where USL contains interactions of all water molecules with the surface and UHB is the system energy rise due to loss of hydrogen bonds at the interface. As discussed above, the bifurcating configuration warrants a lower bound for E(nHB) and thereby UHB can be neglected. Alternatively, USL can be estimated by the repulsive component, since water molecules are mostly distributed above the depletion zone:where ρL is the number density of liquid film. Combining eqn (6) and (8), a quadratic dependence of Wadh on εSL can be derived:

Eqn (9) bridges the gap between the macroscopic interfacial property of Wadh, which is experimentally accessible, and the microscopic parameter εSL of the depth of the Lennard-Jones potential well. Following the standard algorithm to calculate each component of interfacial tensions,[46] we calculated Wadh for different εSL. Fig. 4(b) presents Wadh for different εSL and coincides with our prediction of the quadratic dependence of Wadh on εSL.

Conclusion and perspectives

In this work, we study the interfacial structures of liquid water interacting intimately with a solid surface with varied wettability via MD simulations. Both our simulation results and theoretical analyses unveil the bifurcating configuration of interfacial hydrogen bonds, due to the collective effects of water density depletion, hydrogen bonds and local polarization. Such bifurcating configuration, which is stabilized by the strong confining effect of the solid surface, is comprised of planar water sheets piled up in the perpendicular direction and inter-layer defects of hydrogen bonds bridging two adjacent sheets. As a result, a minimum number of 2.475 hydrogen bonds per water molecule is warranted and a quadratic dependence of work of adhesion on the surface–liquid affinity εSL is disclosed.

In our simulations, the gold surface is essentially set to be electrically neutral and the coulombic interactions between the solid surface and the water film are hence neglected. However, possible surface charge and dielectric polarization may play a certain role in the interfacial properties of water as well, since the dipole moment of water molecules are exceptionally sensitive to the direction of an external electrical field. The introduction of an external electrical field is typically associated with the enhancement of surface wettability, since the electrostatic energy at the solid–liquid interface can be treated as part of the surface free energy. However, it may also alter the water density profile and the formation of hydrogen bonds by tuning the energetically favorable orientation of interfacial mean dipoles. In this respect, the proposed bifurcating configuration can be either stabilized or destabilized, depending on the strength and direction of the external electrical field. The interfacial structures of water molecules, especially the hydrogen bonding formation, on charged and electrically polarizable surfaces or under an external electrical field, warrant further investigation.

Conflicts of interest

There are no conflicts to declare.