Atomistic simulations of friction at an ice − ice interface

Even though the slipperiness of ice is important both technologically and environmentally and often experienced in everyday life, the nanoscale processes determining ice friction are still unclear. We study the friction of a smooth ice–ice interface using atomistic simulations, and especially consider the effects of temperature, load, and sliding velocity. At this scale, frictional behavior is seen to be determined by the lubricating effect of a liquid premelt layer between the sliding ice sheets. In general, increasing temperature or load leads to a thicker lubricating layer and lower friction, while increasing the sliding velocity increases friction due to viscous shear.

One of the many remarkable properties of ice is its low friction coefficient [25][26][27]. It is generally accepted that ice is slippery because of a liquid layer that covers its surface [27][28][29][30], but the mechanism of formation of such a lubricating layer has been a matter of debate and study since the 19th century. In attempts to understand the unusual slipperiness of ice, Faraday (1850) proposed a theory of surface premelting, the spontaneous formation of a liquid layer at the surface of ice well below its melting temperature [28][29][30]. However, the theory soon became controversial-an alternative approach was proposed by Thomson (1857) who formulated the linear dependence between the freezing point depression and applied pressure, and suggested a mechanism of pressure melting as an explanation for liquid layer formation. For many years pressure melting was considered to be the main reason responsible for the low friction coefficient of ice, but later calculations (Bowden and Hughes, 1939) revealed that in standard sliding scenarios the pressure effect is not sufficient to cause surface melting and the biggest contribution comes from frictional heating [31]. However, neither pressure melting (at very low temperatures) nor frictional heating explains why ice can be slippery when one is standing motionless on it.
A number of experimental techniques, such as atomic force microscopy (AFM) [32], nuclear magnetic resonance (NMR) [33,34], X-ray diffraction [35], and photoelectron spectroscopy [36], have been used to study the structural properties of the surface of ice [37]. These experiments provide evidence for the existence of structural disorder at the surface at temperatures below the bulk melting point. The fact that the periodic crystal structure terminates at the surface results in relatively weaker bonding and therefore higher mobility of the surface molecules at temperatures as low as 200 K, consequently the surface molecules show more disordered arrangement [38,39]. However, the temperature range in which this disorder is seen in experiments varies over the techniques applied, as each of them measures different physical properties of the system [37]. Extensive molecular dynamics (MD) simulations have been performed on the surface of ice demonstrating the presence of a quasi-liquid layer at the surface [40][41][42][43][44][45], although it's thickness and temperature dependence is somewhat a function of the potential that is used in simulations to represent the intermolecular interactions.
Despite the importance of ice-ice friction and its relevance to understanding the surface properties of ice, studies of ice friction have been mainly limited to heterogeneous materials [31,[46][47][48][49][50][51][52][53]. The friction of ice on ice has received little attention-the few existing experimental studies show generally low coefficients of friction (from 0.05 up to 1.6) at temperatures close to the melting point, with a clear dependence on the sliding velocity and temperature [54][55][56][57][58][59][60]. In contrast, low temperature studies at high pressure showed minimal dependence on temperature and velocity [61]. Due to the lack of theoretical studies of ice friction, and ice-ice friction in particular, the aim of this work is to systematically study the mechanics of friction between two hexagonal ice (Ih) surfaces using molecular dynamics simulations with a well-established force field. We calculate the frictional properties as a function of temperature, sliding velocity and load, and compare with previous experimental studies. The results can be used as a benchmark for further studies of friction in heterogeneous ice systems and investigations into the role of defects and impurities.

Methods
The system we studied consisted of two parallel slabs of ice Ih, each measuring approximately 3 nm × 3 nm × 3 nm, stacked on top of each other along the z-coordinate of the rectangular simulation box. Both slabs had a similar orientation, exposing the (0001) surface perpendicular to the z-axis. In order to control the distance between the two slabs, and account for the missing macroscopic continuation of the system, the molecules in one of the outer layers in each slab were restrained in a harmonic potential acting in the z direction. Periodic boundary conditions were applied along x, y, and z, and ~10 nm of vacuum was added to the simulation box along z, to minimize spurious interactions with periodic images. In order to simulate friction measurements at constant sliding velocity, the centers of mass of the two slabs were pulled in opposite directions along the x direction, parallel to the surfaces. The harmonic potential in which the centers of mass were pulled had a force constant of 10 4 kJ/(mol•nm 2 ). The harmonic potential used to control the inter-layer distance, as well as the applied load in the friction simulations, was applied to a single layer of water molecules in each slab and had a force constant of 15 × 10 3 kJ/(mol•nm 2 ). All quantities of interest were computed from system configurations saved every 5 ps from 5 ns MD trajectories, following an initial 10 ns MD run to establish surface premelting and equilibration. A snapshot of the system is depicted in Fig. 1.
The two ice slabs were initially well separated in the simulation box and slowly brought in contact by a constant pulling velocity in the z direction (v = 0.001 nm/ps). Once at a predefined distance, measured by the constrained layers, the z coordinates of these layers were fixed and pulling in the x-direction could start. An effective load during sliding was then represented in terms of an average harmonic force acting on the constrained layer perpendicular to the sliding direction and controlled by varying the separation.
The 4-point transferable intermolecular potential (TIP4P) series [62][63][64][65] were used to describe atomistic interactions, and molecular dynamics simulations and data analysis were carried out using the GROMACS (version 4.5.3) simulation package [66]. The Lennard-Jones and short-range electrostatic interactions were truncated at 0.9 nm, and an analytic correction to the dispersion term was applied. The Particle-mesh Ewald (PME) scheme was used to treat the long-range electrostatics. The equations of motion were integrated with the Leapfrog algorithm using a 1 fs timestep. A Nose-Hoover thermostat with a 0.1 ps time constant was applied to the system.
The orientation order parameter proposed by Errington and Debenedetti [67] was used to study the effect of premelting and friction on the ice surfaces. It is a measure of the local tetrahedrality around molecule i, defined as: where ijk Φ is the angle between the oxygen atom of molecule i, and the oxygen atoms on two of its neighbors, j and k. i q takes the value of 1 for a perfect tetrahedral structure and 0 for complete disorder. In practice, the order parameter is about 0.95 in bulk ice and of the order 0.5−0.85 in liquid water.

Results
To first establish the surface premelting behavior in simulations, we studied the nature of a free (0001) surface. For comparison, we used the water models simple point charge (SPC), TIP4P, TIP4P/2005, and TIP4P/Ice with melting temperatures of 190 K, 232.0 K, 252.1 K and 272.2 K, respectively [44]. Due to the differences in melting points, also the premelting temperatures (temperature at which the disordered layer thickness reaches 1 Å [44]) were seen to depend on the potential, but in all cases the thickness of  [36,44]. Figure 2 shows the order parameter, calculated as a function of surface depth, for the SPC and TIP4P/Ice water potentials at temperatures below their melting points illustrating the similarity in premelting behavior. The number density profiles calculated for the TIP4P/Ice system at minimum (230 K) and maximum (270 K) temperatures indicate a layered structure of the premelt (see Fig. 3). In addition, density profiles calculated separately for the bulk and the surface layers at 230 K and 270 K temperatures show the difference between the crystalline order and the surface disorder (see Fig. 4), more pronounced close to the melting point (see Figs. 4(b) and 4(d)). Due to the different packing orders of the crystal, in the direction perpendicular and parallel to the sliding interface (ice−vacuum interface), the calculated average densities in Figs. 3 and 4 differ significantly. We chose TIP4P/Ice, which matches the experimental melting temperature [64], for further   studies of friction and do not expect the results to depend significantly on the flavor of the TIP4P model beyond the temperature scale.
When the formation and growth of the liquid layer at the surface of ice is greatly influenced by the temperature then its frictional behavior is also expected to be temperature dependent. Figure 5(a) shows the temperature dependence of frictional force for a "moderate load" case. The lowering of friction as temperature rises can be understood in terms of increased lubrication. Generally, the molecules at the interface will form hydrogen bonds between each other resisting sliding, but in the premelt the molecules are less coordinated (q is lower, see Fig. 5(b)) and so are also more weakly bound. The thickness of the premelt layer increases with temperature resulting in better lubrication at the ice−ice interface. Also, due to the increased thermal motion of the molecules at higher temperatures the average hydrogen bond strength effectively weakens [68].
Similarly also the sliding velocity influences friction, since it is directly related to the frictional heating and therefore plays significant role in the resulting lubrication. At higher sliding velocities more frictional heat is generated, increasing the thickness of the interfacial water layer. During sliding, frictional heating will locally raise the temperature at the contact layer. However, in order to separate the effects of heat and sliding velocity, the temperature of the liquid layer was also kept approximately constant with the thermostat during simulations. Friction is seen to increase linearly with increasing sliding velocity, which can be interpreted to be due to viscous shear in the liquid layer between sliding surfaces (see Fig. 6).
Next, we examine the effect of load on frictional force. As the two ice slabs are brought together, the thickness of the liquid layer between them decreases initially. This happens because when the ice sheets are far apart, in the negative load regime, water molecules fill the small void between the surfaces. As the slabs are pressed together, the liquid layer is at first confined in a smaller space, increasing the density, until at a high load the amount of liquid starts to increase due to pressure melting. At this point the thickness of the liquid layer starts to increase again. Furthermore, the diffusion constant (computed using the Einstein relation) of water molecules confined at the interface decreases with the slab separation (see Fig. 7). The diffusion constants for each temperature and separation distance were calculated at sliding conditions in order to avoid freezing of the interface, especially for lighter loads. The correspondence between effective load and slab separation is shown in Fig. 8. Also here it is apparent that we cross from attractive regime (negative load) to repulsive regime (positive load) when the constraining force changes sign. The effect of load on frictional force is shown in Figs. 9 and 10 in tandem    with both the sliding velocity (see Fig. 9) and temperature (see Fig. 10). In all cases, in the attractive regime, frictional force decreases weakly as the slabs are brought closer. Near the crossover to repulsive regime, we see a rapid drop of frictional force followed again by fairly weak load dependence in the repulsive regime. Finally, by transforming from separation to load using the dependence of Fig. 8, we can calculate the coefficient of friction. This is shown in Fig. 11 where the friction coefficient is seen to decrease linearly with temperature and increase with velocity. The coefficient is also found to decrease as the applied load increases (see Fig. 12), whereas the effect of adhesion is seen at lighter loads.    12 Coefficient of friction as a function of separation distance between harmonically restrained layers (applied load) at 240 K, 250 K, and 260 K temperatures. The coefficient increases rapidly around 5.5-5.6 nm separation where the load goes from attractive to repulsive regime as the ice slabs are pressed together. This can be explained by the attraction between interface water molecules at smaller loads resisting sliding.

Discussion
Experimental studies have shown that generally friction on ice surfaces is influenced by temperature, sliding velocity and applied load [31,[46][47][48]54]. However, when comparing the presented simulations to experiments, it should be noted that the simulated interfaces are atomistically smooth and always separated by the liquid premelt layer. This suggests we are simulating a hydrodynamic friction regime, where friction is due to viscous shear of the liquid film, corresponding to experiments done close to the melting point where thick premelt layers are expected.
Pure viscous shearing would imply a frictional force proportional to the sliding velocity. However, the simulations show there is also a temperature dependent static friction component-a minimum force needed to move the ice sheets at low velocities, making the premelt act more like a Bingham plastic [69] rather than a Newtonian liquid. Simulations were also carried out with a constant sliding force, chosen to be lower than the minimum force calculated from constant velocity simulations, verifying that a finite force is necessary to initiate sliding. This is explained by the tendency of the interface liquid layer to solidify between the ice sheets. The lower the temperature is, the stronger the pulling force needs to be to prevent the system from freezing to a single piece of solid ice.
Experimentally, friction between ice surfaces decreases with increasing sliding velocity due to increased frictional heating and thicker liquid layer at the interface. However, close to the melting point, the coefficient of friction becomes proportional to v 1/2 [54] once the interface is completely covered in liquid and viscous shear becomes dominant. For temperatures close to the melting point, our simulations show a similar dependence (The relative root mean square error and the correlation coefficient of fitting is 0.07 and 0.1 respectively) of the frictional force on sliding velocity. At low temperatures the difference between experimental and simulation results can be understood to be due to the surface roughness present in experimental systems, missing from the simulations.
Ice friction experiments performed at the macroscopic scale show a decrease in the friction coefficient with increasing normal load, with a considerable difference between results depending on the material sliding over ice surface as well as temperature and velocity [31,47,48,52,54]. Similar trends are also seen in the simulations, which also show a decrease in the frictional force with increasing load, although in all cases (regardless of temperature or sliding velocity) the dependence becomes less pronounced as the load increases.