Dynamics of nanoscale droplets on moving surfaces.

We use molecular dynamics (MD) simulations to investigate the dynamic wetting of nanoscale water droplets on moving surfaces. The density and hydrogen bonding profiles along the direction normal to the surface are reported, and the width of the water depletion layer is evaluated first for droplets on three different static surfaces: silicon, graphite, and a fictitious superhydrophobic surface. The advancing and receding contact angles, and contact angle hysteresis, are then measured as a function of capillary number on smooth moving silicon and graphite surfaces. Our results for the silicon surface show that molecular displacements at the contact line are influenced greatly by interactions with the solid surface and partly by viscous dissipation effects induced through the movement of the surface. For the graphite surface, however, both the advancing and receding contact angles values are close to the static contact angle value and are independent of the capillary number; i.e., viscous dissipation effects are negligible. This finding is in contrast with the wetting dynamics of macroscale water droplets, which show significant dependence on the capillary number.


■ INTRODUCTION
Wetting phenomena play an important role in diverse processes across physics, chemistry, and biology. 1−3 The physics of wetting phenomena for water droplets is of fundamental importance in the design of surfaces that can mimic natural surfaces. 4−6 A thorough understanding of solid−liquid interactions at a molecular level is also crucial to technological applications, including surface coating, emulsions, oil recovery, and in microfluidic and nanofluidic applications. 7−10 The key experimental parameter describing the degree of wetting is the static contact angle θ s , measured through the liquid L placed in contact with a solid S, at the contact line. The extent to which a liquid wets a given solid, i.e., the wettability, has been a subject of much research over the past few decades, both theoretically and experimentally. 11 The wettability determines the equilibrium configuration of the system: if θ s is zero, then the liquid is said to wet the solid completely and the solid surface is fully hydrophilic; if it is 180°, then the system is said to nonwet the solid and the surface is fully hydrophobic.
During the dynamic wetting process, in which the contact line moves across a solid surface, there may be several different modes of dissipation. For partially wetting and nonwetting Newtonian liquids (and leaving aside dissipations within the main bulk of the liquid) the principal losses are viscous ones due to flow very near the contact line and contact-line friction associated with the creation (or loss) of the solid−liquid interface. 12 The actual physics that governs the wetting dynamics at the liquid−solid interface remains poorly understood. 13,14 One reason for this is that the dynamics are dictated by physical phenomena taking place on different length scales. The large-scale dynamics are typically governed by hydro-dynamic theory, while the movement of the contact line itself is determined by processes on (or just above) molecular length scales. Another reason why dynamic wetting has remained unclear is that experiments are quite difficult, with a large span in length scales and very rapid time scales. 15 There exist three principal theories for the description of dynamic wetting phenomena, namely, molecular-kinetic theory, 13,12 hydrodynamic theory, 16,17 and phase field theory. 18−21 Molecular-kinetic theory describes dynamic wetting as the disturbance of adsorption equilibria at the contact line. The movement of the contact line is determined by the statistical dynamics of the molecular motion within the three phase zone, where the solid, liquid, and gas phases meet. Hydrodynamic theories are built on a continuum description and typically use the lubrication approximation involving either a microscopic cutoff length beyond which the solution is truncated or a postulated region of local slip between the liquid and solid. A different way to handle moving contact lines without violating the no-slip boundary condition is the phase field theory, which enables the contact line motion by diffusive interfacial fluxes. Its theoretical framework stems from a thermodynamic formulation, 21 based on a description of the free energy in the system. Although all three theoretical formulations have been applied with some success, the models involve adjustable parameters that need to be determined through fitting experimental data. In addition, the applicability of these models to nanoscale droplet dynamics is questionable, as the exact nature of hydrophobic and hydrophilic interactions at the nanoscale remains elusive and controversial. 2,7−9 Molecular dynamics (MD) simulations have been extensively used to study the surface wetting phenomena of nanoscale droplets. 22−28 MD is found to be an accurate deterministic method based on classical equations of motion and allowing for realistic molecular behavior, i.e., molecular attractions, repulsions, movements, and scatterings. Microscopic equivalence of the contact angle has then often been casually applied to nanometer-scale droplets on a wide variety of surfaces, such as polymers, cellulose, silica, graphite, and carbon nanotubes, among others, in order to establish a connection between the microscopic calculations and the macroscopic wetting properties of the surfaces. The equilibrium contact angle has been used as a reference in order to tune the intermolecular interaction parameters. In most of the MD-based studies, the surfaces have been held frozen during these simulations in order to save computational time.
Sergi et al. 29 performed MD simulations of water droplets on graphite surfaces with an alternative approach to measuring the contact angle. Wang and Zhao 30 investigated the contact angle hysteresis of nanodroplets on both rigid and flexible substrates with different wettabilities, in response to a body force. Koishi et al. 31 and Savoy et al. 32 carried out simulations of the wetting characteristics of water and oil droplets, respectively, on static pillar-type nanostructures, and Jeong et al. 33 performed similar studies by taking into account the effect of a body force on the water droplet behavior. Dutta et al. 34 studied size and temperature effects on the wetting transition of water on graphite and boron nitride surfaces.
Previous MD investigations have examined the effect on the wetting behavior of temperature, size, surface physical chemistry, gravity-like forces, and nanostructure effects. However, none of the MD investigations considered wall movement effects on the wetting dynamics of water droplets. However, these conditions are relevant to many industrial applications, such as slot bead, slot curtain, and roller coating methods, 10,35−37 and other experimental studies. 38 In this paper we report investigations on the interactions of water droplets with moving silicon and graphite surfaces at a nanometer length scale. We use the molecular dynamics (MD) simulation method. Two different surfaces are studied, and their velocities are varied in order to gain insights into the correlations between the surface movement and the contact angle. Advancing and receding contact angles, and contact angle hysteresis measurements, are obtained at different capillary numbers for silicon and graphite surfaces. The structure of the nanoscale droplets is also studied by examining the water density profiles, water depletion layers, and hydrogen bonding near both the hydrophobic (graphite) and hydrophilic (silicon) surfaces. Water droplet behavior on a fictitious superhydrophobic surface is also investigated under static conditions. The existence of a vapor-like depletion layer of water molecules near a superhydrophobic surface is also discussed.

■ METHODOLOGY
To study the dynamic wetting of water on moving surfaces, MD simulations of water droplets on smooth graphite and silicon are performed. Silicon is a hydrophilic material while graphite is weakly hydrophobic. In the following, the MD technique is described, along with details of how the dynamic contact angle, contact angle hysteresis, density profiles, and hydrogen bond distribution are measured from the numerical experiments.
Molecular Dynamics Simulations. We use the rigid TIP4P/2005 water model 39 to model the condensed phases of water. Studies have shown that this model can reproduce major water properties more accurately than other commonly used rigid models, such as the SPC/E, TIP3P, and TIP4P. 40 It consists of a Lennard-Jones (LJ) center on the oxygen (O) with ε OO = 0.7749 kJ mol −1 and σ OO = 0.315 89 nm and three fixed point charges on the fourth massless M site (−1.1128 e) and the two hydrogen atoms (0.5564 e). Force calculations due to Coulomb interactions are truncated at the same cutoff distance of 1.2 nm as the LJ interactions. The electrostatic interactions are also shifted, a common practice in MD simulations, which offers adequate accuracy and significant computational savings compared to other methods. 26,41−43 In this study, Hamilton's quaternions 44 where U s is the solid wall velocity, μ = 0.855 × 10 −3 Pa/s is the dynamic macroscopic viscosity of water, and γ = 69.3 × 10 −3 N/m is the surface tension according to the detailed study of Vega and Abascal. 40 No forces are applied to the graphite or silicon atoms, but their relative distance with each other is kept fixed to represent an inert wall; solid atoms are only allowed to move as a bulk in the direction of the given velocity. It has been shown previously that fixing the solid atoms of a surface does not affect the contact angle of the droplets, but it does reduce significantly the computational expense. 26 The dimensions of the surfaces are 20 × 20 × 0.34 nm and 38.018 × 38.018 × 1.5 nm for graphite and silicon, respectively. The graphite consists of two staggered graphene sheets with an interlayer distance of 0.34 nm, while the silicon wall is constructed as a uniform crystallite. This silicon model is later adapted to represent a fictitious hydrophobic material by altering the interaction parameter only. The thicknesses of the solid surfaces are intentionally kept small, as additional layers of solid material are not expected to have a significant influence on the water due to the employed cutoff radius of 1.2 nm. Periodic boundary conditions are applied in all directions of the simulation box, which in effect means that the droplet lies on a surface of infinite extent. In all our simulations the water−surface interaction is solely based on a Lennard-Jones potential between the oxygen atom of the water molecule and the carbon or silicon atoms of the surface. The LJ values of σ CO = 0.319 nm and σ SiO = 0.323 nm are employed for all our simulations, following the Lorentz− Berthelot mixing rules and previous studies. 26,45,46 The LJ interaction parameters ε CO = 0.427 kJ mol −1 , ε SiO = 2.36 kJ mol −1 , and ε fiO = 0.139 kJ mol −1 have been chosen in order to reproduce the macroscopic static contact angle of the simulated water model on the selected surfaces, namely 86°for graphite, 47 43°for silicon, 48 and 150°for the fictitious hydrophobic material.
Our simulations are performed using mdFoam, 49 simulations are carried out for a problem time of 1 ns, with an integration step of 2 fs. Contact angle, density, and hydrogenbonding measurements are averaged over a minimum of 500 000 samples, which is sufficient to obtain accurate results within ±1% error. The simulations are run for 200 ps in order for the molecular ensemble to relax before sampling begins. During this relaxation time the system is coupled to a Berendsen thermostat at a temperature of 300 K, which is then removed for the remainder of the simulation and subsequent measurements. The sampling is performed in the microcanonical ensemble (NVE) of constant number of atoms, constant volume, and constant energy. The evaporation time of the smallest droplet we simulate is of the order of 100 μs, 52 which is quite large compared to the simulation time scale.
The water molecules are placed initially in a rectangular lattice of dimensions 5.8 × 6.2 × 1.6 nm that contains 2000 water molecules for the graphite surface or a 6.3 × 3.8 × 6.3 nm rectangular lattice that contains 5000 water molecules for the silicon surface. This difference is necessary because the hydrophilic nature of the silicon surface means that 2000 water molecules for a droplet are insufficient to provide the necessary thickness for accurate contact angle measurement, as explained in the following section. Samples of the molecular trajectories and other measurements are stored every 0.2 ps. Our simulations typically took 5 days each, using 4 or 8 cores of a parallel computer for the graphite and silicon surfaces, respectively.
Measurement Techniques. Water density isochore profiles are obtained from the MD simulation trajectories. To do this, the volume around the water droplet is divided into bins radially from the center of mass of the droplet and also in the direction perpendicular to the solid surface. This results in a cylindrical binning, with the center of mass of the water droplet as the reference axis' origin. The bins have equal volume, and the average bin height is 0.05 nm.
In order to capture the dynamic contact angle, the cylindrical binning follows the center of mass of the water droplet at every time step. In each bin the water density is calculated, and in addition, the binning cylinder is divided into two parts: one for the advancing part of the droplet and one for the receding part (depending on the direction of the imposed velocity of the wall). With the resulting contour plots, a two-step procedure is followed in order to calculate the advancing (θ) and receding (ϕ) contact angles, similar to that described in ref 53. First, the liquid−vapor interface contour line is chosen as the one with half the density of bulk water. Second, a circular best fiti s applied on these points and is extrapolated to the solid surface, where the contact angles θ and ϕ are then measured. It should be noted that density contour points at a distance less than 0.8 nm from the solid surface are not taken into account during the fitting. This ensures that any fluctuations at the liquid−solid interface do not affect the measurements. It was shown in previous studies that the choice of this distance has negligible effect on the measured contact angles. 26 Hence, we decided to perform droplet simulations on silicon with 5000 water molecules. On the graphite sheet, water does not spread so much, which leads to a thick droplet formation, even with 2000 molecules. The difference between the two measured angles then represents the contact angle hysteresis (H = θ − ϕ). All density values, lengths, and angles presented in the following figures are in reduced units, with a reference density of 952.03 kg/m 3 , a reference length of 0.3154 nm, and the static contact angle of each material (θ s ).
The hydrogen bonds (HB) distribution in the droplet is calculated using a similar binning procedure as for the water density. The geometrical criterion of Luzar and Chandler 54 is used to count the HB: two water molecules are considered to be hydrogen bonded if the oxygen and the hydrogen that form the HB are less than 0.35 nm apart and if the oxygen−oxygen/ oxygen−hydrogen bond angle is less than 30°.
Size effects are important when contact angle measurement is considered. The modified Young's equation 55 correlates the microscopic contact angles θ and ϕ (or θ D when we refer to a general microscopic contact angle) to the macroscopic ones, θ ∞ and ϕ ∞ , respectively. In the following figures and results, the contact angles are always the microscopic ones.

■ RESULTS AND DISCUSSION
We first make qualitative observations on the droplet behavior before discussing the quantitative results. Figures 1 and 2 show typical snapshots and density contours, respectively, of water droplets at their steady state on three different static surfaces: (a) silicon, (b) graphite sheet, and (c) the artificial hydrophobic surface. In Figure 1, the water molecules comprise red (oxygen) and white (hydrogen) atoms; the solid surface molecules are in brown (silicon), gray (graphite), and green (artificial hydro- phobic surface). The water droplet relaxes from its initial rectangular configuration to roughly a capped spherical shape within the first 100 ps in each simulation and remains so for the rest of the simulation. The shape of the droplet at equilibrium is almost flat on the silicon surface and hemispherical on the graphite. The water droplet adopts a near-perfect spherical shape on the artificial hydrophobic surface and completely nonwets the surface from the molecular point of view. We also observe that the position of the water droplet remains almost unchanged during the simulation with the silicon surface, while the droplet becomes more mobile laterally as the hydrophobicity increases. From the density contours in Figure 2, static microscopic contact angles for the silicon, graphite, and the artificial hydrophobic surfaces are θ s = 26.8°, 88.8°, and 180°, respectively.
In Figure 3, we show the density profiles of water droplets on the three different surfaces along the direction normal to the surface (i.e., the z-direction). The nanoscale droplet completely spreads on the silicon surface, which results in local oscillations in the density profile without reaching a constant value in the bulk. Oscillations around the bulk value for the first two or three hydration shells from the surface reflect a high degree of spatial ordering of water molecules in this region. A density higher than the bulk value for z < 0.95 also implies that the first water layer may form a relatively tight boundary at the interface. The water density profile on the graphite surface follows the same trend as the silicon case for z < 0.95. At larger distances from the surface, the water density levels off to the uniform bulk density, which is a typical indication of a random isotropic distribution of the molecules. For the case of the strongly hydrophobic surface, water molecules are not seen at all up to z = 9.5, and for higher values of z the density increases monotonically to reach its bulk value.
The water depletion layer, which we define as the region between the surface and the height where the water density falls below half the bulk density, can also be evaluated from Figure  3. The layer thickness is negligible, less than 0.95 for both the silicon and graphite surfaces (around the size of two water molecules). This is consistent with the results of other recent studies showing that the thickness of the depletion zone can be a few angstroms for both hydrophilic and weakly hydrophobic surfaces. 56 However, for the superhydrophobic surface, the thickness of the water depletion layer is found to be around 12, which is well outside the interaction range of water−solid molecules. This depletion layer causes liquid droplets to effectively slip over hydrophobic surfaces and is observed experimentally. 57 In general, the existence of a vapor layer near the interface is more likely for a hydrophobic surface. This is primarily because the structure of water molecules next to a hydrophobic surface is less ordered than in the bulk phase, while the cohesive strength of water may be significantly reduced. This is also consistent with Figure 4, where the number of hydrogen bonds per water molecule close to the hydrophobic surface is shown to be zero. Once the nucleation barrier is exceeded, either the incipient depletion layer formed on the solid surface or, in more realistic conditions, gas molecules may be trapped in the gap between the liquid and solid layers. This effect acts to prevent the liquid from being directly exposed to the wall surface. In such cases, the liquid is not likely to experience the presence of the wall directly and may smoothly "sail over" the intervening depletion layers, instead of being in proximate contact with the wall. Figure 4 shows the average number of hydrogen bonds per water molecule n HB along the direction normal to the surface for the different surfaces. The profiles for silicon and graphite surfaces reveal that the first peak in n HB is located between the  first and the second peaks in the water density profile. The number of hydrogen bonds drops quickly near the surface as the very first layer of water molecules interfacing with the surface cannot form hydrogen bonds in the direction toward the surface. The rise in the number of hydrogen bonds between the first and second layer of water molecules indicates a locally tight structure of these two water layers. For the artificial hydrophobic surface case, the profile increases monotonically before reaching its constant bulk value. Figure 5 presents the variation of the normalized advancing and receding dynamic contact angles with capillary number (Ca) for water droplets on moving silicon and graphite surfaces. The majority of the measured angles on the moving silicon surface (Figure 5a) are smaller than the static one (θ s ). For Ca rising between 0.01 and 0.1, the advancing contact angle reduces from 1.04 to 0.78. It then increases to 0.93 by Ca = 0.2 and then decreases slowly with further increases in Ca. The receding contact angle increases from 0.88 to 1.02 as Ca increases from 0.01 to 0.025, but for higher values of Ca it fluctuates between 0.85 and 0.95. On the other hand, Figure 5b indicates that the contact angle of nanoscale droplets on graphite surfaces does not show a large variation with the capillary number. Independent of the velocity of the graphite wall, both the dynamic advanced and receding contact angles are close to the static contact angle value. Figure 6 shows the normalized contact angle hysteresis variation with Ca for water droplets on moving surfaces. Normalized hysteresis values for the graphite surface are close to 0°, with only small fluctuations. For the silicon surface, however, hysteresis values decrease from 0.17 to −0.07 as Ca increases from 0.01 to 0.1, and they increase to 0.02 with further increase in Ca.
For the silicon surface, Figures 5 and 6 indicate that we can classify the Ca range into three distinctive regimes: the first for Ca smaller than 0.1, the second for values between 0.1 and 0.2, and a third for higher Ca values. In the first regime, surface tension forces are dominated by viscous stresses that are induced through the boundary movements. In the second regime, both viscous and inertial forces are equally important, and in the third regime inertial forces dominate the other two effects. For the graphite surface, molecular displacements at the contact line are purely influenced by interactions with the solid surface (i.e., surface tension forces), and both the viscous dissipation effects and inertial forces are negligible. Figures 7 and 8 show the normalized water density profiles with normalized distance from the surfaces (z-direction) for silicon and graphite, respectively. In the case of silicon, the profiles are nearly identical for all capillary numbers. However, deviations can be noticed around the first peak in the density profiles. At a given z, for the range 0.01 < Ca < 0.1, density increases with an increase in Ca, while it decreases for Ca > 0.1. Denser packing of water molecules in the first layer causes a   relatively tight boundary at the interface, which may lead to a decrease in the contact angle. This finding is also consistent with the measurements in Figure 5a, where the contact angle decreases from Ca = 0.01 to 0.1 but increases with further increase in Ca. In the case of graphite, the density profiles are negligibly affected by variations in the capillary number, although minor deviations can be noticed around the first peak and crest. Figures 9 and 10 show the average number of hydrogen bonds per water molecule along the direction normal to the silicon and graphite surfaces, respectively. For silicon, there are no hydrogen bonds up to a normalized distance of 0.35 from the surface. For Ca = 0.01, the number of hydrogen bonds then increases monotonically to its bulk value. For the Ca = 0.025, 0.05, and 0.2 cases, however, a rise in the number of hydrogen bonds is seen before the first layer (z < 0.8) of water molecules. This means a relatively tighter binding of water molecules with the surface through the extra hydrogen bonds in this near-wall region. This may lead to better wetting conditions for water droplets and a decrease in the contact angle. A further increase in Ca leads again to a monotonic increase in the number of hydrogen bonds to its bulk value, without any local peak near the surface. In the case of graphite, the number of hydrogen bonds is zero up to a normalized distance of 0.7, i.e., twice that for the silicon surface. This means there are poor wetting conditions near the surface, irrespective of the capillary number. Viscous forces have a negligible effect on the number of   hydrogen bonds in both the first and second layers (i.e., 0.8 < z < 2) of water molecules and in the bulk region (z > 2).

■ CONCLUSIONS
We have performed molecular dynamics (MD) simulations of nanoscale water droplets on different static and moving surfaces. Silicon, graphite, and an artificial material modeled hydrophilic, weakly, and strongly hydrophobic surfaces, respectively. Water density and hydrogen-bonding profiles were extracted, and both the effect of solid−liquid interactions on the formation of water depletion layers and the cohesive strength of water molecules have been discussed. A large depletion layer is formed near the hydrophobic surface, primarily because the structure of water molecules next to a hydrophobic surface is less ordered than in the bulk phase, which significantly reduces the cohesive strength of the water.
In addition to the solid−liquid interaction effects, viscous dissipation effects were also investigated by moving the silicon and graphite surfaces. It was found that for nanoscale droplets the solid−liquid interactions play a vital role in determining the wetting dynamics, while viscous dissipation effects induced by the moving surface were found to be only slightly important for the silicon surface and negligible for the graphite surface. These observations are different from the wetting dynamics of macroscale droplets, which show a significant dependence on the capillary number. For a silicon surface, the advancing contact angle decreases with an increase in capillary number through 0.01 < Ca < 0.1 and increases for 0.1 < Ca < 1. In the former regime, surface tension plays a key role, while in the latter both viscous and inertial forces are important. In the case of the graphite surface, solid−liquid interactions play a major role irrespective of capillary number.
Our study has yielded physical insight into the near-surface wetting dynamics of nanoscale droplets, which may aid in developing new boundary conditions for continuum and hybrid models of liquid flows in micro-and nanoscale devices. For example, it is possible to apply molecular kinetic theory to our MD data in order to extract the coefficient of contact line friction, the molecular jump frequency, and the molecular jump length. The effect of droplet size on the measured properties can then be reported as a function of liquid Knudsen number, 20 which is the ratio of the molecular jump length to the characteristic macroscopic length. Other future work should include (i) obtaining measurements for more capillary numbers between 0.01 and 1 for various different surfaces and (ii) repeating the simulations for larger droplets in the whole range of capillary numbers in order to investigate size effects (although in this case the MD simulations will be exceptionally time-consuming).

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
The research leading to these results has received funding from the UK's Engineering and Physical Sciences Research Council under Grant EP/I011927/1. Results were obtained using the ARCHIE-WeSt High Performance Computer, under EPSRC grant EP/K000586/1.