ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceWater2/iceWater3.tex
Revision: 4255
Committed: Mon Dec 15 17:10:05 2014 UTC (10 years, 7 months ago) by plouden
Content type: application/x-tex
File size: 50219 byte(s)
Log Message:
Added the SPC/E Eta(225K) shear rates.

File Contents

# User Rev Content
1 gezelter 4243 %% PNAStwoS.tex
2     %% Sample file to use for PNAS articles prepared in LaTeX
3     %% For two column PNAS articles
4     %% Version1: Apr 15, 2008
5     %% Version2: Oct 04, 2013
6    
7     %% BASIC CLASS FILE
8     \documentclass{pnastwo}
9    
10     %% ADDITIONAL OPTIONAL STYLE FILES Font specification
11    
12     %\usepackage{PNASTWOF}
13     \usepackage[version=3]{mhchem}
14 gezelter 4245 \usepackage[round,numbers,sort&compress]{natbib}
15     \usepackage{fixltx2e}
16 gezelter 4247 \usepackage{booktabs}
17     \usepackage{multirow}
18 gezelter 4245 \bibpunct{(}{)}{,}{n}{,}{,}
19     \bibliographystyle{pnas2011}
20 gezelter 4243
21     %% OPTIONAL MACRO DEFINITIONS
22     \def\s{\sigma}
23     %%%%%%%%%%%%
24     %% For PNAS Only:
25     %\url{www.pnas.org/cgi/doi/10.1073/pnas.0709640104}
26     \copyrightyear{2014}
27     \issuedate{Issue Date}
28     \volume{Volume}
29     \issuenumber{Issue Number}
30     %\setcounter{page}{2687} %Set page number here if desired
31     %%%%%%%%%%%%
32    
33     \begin{document}
34    
35 gezelter 4254 \title{The different facets of ice have different hydrophilicities:
36     Friction at water / ice-I\textsubscript{h} interfaces}
37 gezelter 4243
38     \author{Patrick B. Louden\affil{1}{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame,
39     IN 46556}
40     \and
41     J. Daniel Gezelter\affil{1}{}}
42    
43     \contributor{Submitted to Proceedings of the National Academy of Sciences
44     of the United States of America}
45    
46     %%%Newly updated.
47     %%% If significance statement need, then can use the below command otherwise just delete it.
48     %\significancetext{RJSM and ACAC developed the concept of the study. RJSM conducted the analysis, data interpretation and drafted the manuscript. AGB contributed to the development of the statistical methods, data interpretation and drafting of the manuscript.}
49    
50     \maketitle
51    
52     \begin{article}
53 gezelter 4245 \begin{abstract}
54 gezelter 4254 We present evidence that some of the crystal facets of
55     ice-I$_\mathrm{h}$ posess structural features that can reduce the
56     effective hydrophilicity of the ice/water interface. The spreading
57     dynamics of liquid water droplets on ice facets exhibits long-time
58     behavior that differs substantially for the prismatic
59 gezelter 4245 $\{1~0~\bar{1}~0\}$ and secondary prism $\{1~1~\bar{2}~0\}$ facets
60     when compared with the basal $\{0001\}$ and pyramidal
61     $\{2~0~\bar{2}~1\}$ facets. We also present the results of
62     simulations of solid-liquid friction of the same four crystal
63     facets being drawn through liquid water. Both simulation
64     techniques provide evidence that the two prismatic faces have an
65     effective surface area in contact with the liquid water of
66     approximately half of the total surface area of the crystal. The
67     ice / water interfacial widths for all four crystal facets are
68     similar (using both structural and dynamic measures), and were
69     found to be independent of the shear rate. Additionally,
70     decomposition of orientational time correlation functions show
71     position-dependence for the short- and longer-time decay
72     components close to the interface.
73 gezelter 4243 \end{abstract}
74    
75     \keywords{ice | water | interfaces | hydrophobicity}
76     \abbreviations{QLL, quasi liquid layer; MD, molecular dynamics; RNEMD,
77     reverse non-equilibrium molecular dynamics}
78    
79 gezelter 4245 \dropcap{S}urfaces can be characterized as hydrophobic or hydrophilic
80     based on the strength of the interactions with water. Hydrophobic
81     surfaces do not have strong enough interactions with water to overcome
82     the internal attraction between molecules in the liquid phase, and the
83     degree of hydrophilicity of a surface can be described by the extent a
84 gezelter 4254 droplet can spread out over the surface. The contact angle, $\theta$,
85     formed between the solid and the liquid depends on the free energies
86     of the three interfaces involved, and is given by Young's
87     equation~\cite{Young05},
88 gezelter 4245 \begin{equation}\label{young}
89     \cos\theta = (\gamma_{sv} - \gamma_{sl})/\gamma_{lv} .
90     \end{equation}
91     Here $\gamma_{sv}$, $\gamma_{sl}$, and $\gamma_{lv}$ are the free
92 gezelter 4254 energies of the solid/vapor, solid/liquid, and liquid/vapor interfaces,
93 gezelter 4245 respectively. Large contact angles, $\theta > 90^{\circ}$, correspond
94     to hydrophobic surfaces with low wettability, while small contact
95     angles, $\theta < 90^{\circ}$, correspond to hydrophilic surfaces.
96     Experimentally, measurements of the contact angle of sessile drops is
97     often used to quantify the extent of wetting on surfaces with
98     thermally selective wetting
99 gezelter 4254 characteristics~\cite{Tadanaga00,Liu04,Sun04}.
100 gezelter 4243
101 gezelter 4245 Nanometer-scale structural features of a solid surface can influence
102     the hydrophilicity to a surprising degree. Small changes in the
103     heights and widths of nano-pillars can change a surface from
104     superhydrophobic, $\theta \ge 150^{\circ}$, to hydrophilic, $\theta
105 plouden 4246 \sim 0^{\circ}$.\cite{Koishi09} This is often referred to as the
106 gezelter 4245 Cassie-Baxter to Wenzel transition. Nano-pillared surfaces with
107     electrically tunable Cassie-Baxter and Wenzel states have also been
108 gezelter 4254 observed~\cite{Herbertson06,Dhindsa06,Verplanck07,Ahuja08,Manukyan11}.
109 gezelter 4245 Luzar and coworkers have modeled these transitions on nano-patterned
110 gezelter 4254 surfaces~\cite{Daub07,Daub10,Daub11,Ritchie12}, and have found the
111 gezelter 4245 change in contact angle is due to the field-induced perturbation of
112 gezelter 4254 hydrogen bonding at the liquid/vapor interface~\cite{Daub07}.
113 gezelter 4245
114     One would expect the interfaces of ice to be highly hydrophilic (and
115     possibly the most hydrophilic of all solid surfaces). In this paper we
116     present evidence that some of the crystal facets of ice-I$_\mathrm{h}$
117 gezelter 4254 have structural features that can reduce the effective hydrophilicity.
118 gezelter 4245 Our evidence for this comes from molecular dynamics (MD) simulations
119     of the spreading dynamics of liquid droplets on these facets, as well
120     as reverse non-equilibrium molecular dynamics (RNEMD) simulations of
121     solid-liquid friction.
122    
123     Quiescent ice-I$_\mathrm{h}$/water interfaces have been studied
124     extensively using computer simulations. Haymet \textit{et al.}
125     characterized and measured the width of these interfaces for the
126     SPC~\cite{Karim90}, SPC/E~\cite{Gay02,Bryk02},
127     CF1~\cite{Hayward01,Hayward02} and TIP4P~\cite{Karim88} models, in
128     both neat water and with solvated
129     ions~\cite{Bryk04,Smith05,Wilson08,Wilson10}. Nada and Furukawa have
130     studied the width of basal/water and prismatic/water
131     interfaces~\cite{Nada95} as well as crystal restructuring at
132     temperatures approaching the melting point~\cite{Nada00}.
133    
134 gezelter 4243 The surface of ice exhibits a premelting layer, often called a
135 gezelter 4245 quasi-liquid layer (QLL), at temperatures near the melting point. MD
136     simulations of the facets of ice-I$_\mathrm{h}$ exposed to vacuum have
137     found QLL widths of approximately 10 \AA\ at 3 K below the melting
138 gezelter 4254 point~\cite{Conde08}. Similarly, Limmer and Chandler have used the mW
139 gezelter 4245 water model~\cite{Molinero09} and statistical field theory to estimate
140 gezelter 4254 QLL widths at similar temperatures to be about 3 nm~\cite{Limmer14}.
141 gezelter 4243
142 gezelter 4245 Recently, Sazaki and Furukawa have developed a technique using laser
143     confocal microscopy combined with differential interference contrast
144     microscopy that has sufficient spatial and temporal resolution to
145     visulaize and quantitatively analyze QLLs on ice crystals at
146 gezelter 4254 temperatures near melting~\cite{Sazaki10}. They have found the width of
147 gezelter 4245 the QLLs perpindicular to the surface at -2.2$^{o}$C to be 3-4 \AA\
148     wide. They have also seen the formation of two immiscible QLLs, which
149 gezelter 4254 displayed different dynamics on the crystal surface~\cite{Sazaki12}.
150 gezelter 4243
151 gezelter 4245 There is now significant interest in the \textit{tribological}
152     properties of ice/ice and ice/water interfaces in the geophysics
153     community. Understanding the dynamics of solid-solid shearing that is
154 gezelter 4254 mediated by a liquid layer~\cite{Cuffey99, Bell08} will aid in
155 gezelter 4245 understanding the macroscopic motion of large ice
156 gezelter 4254 masses~\cite{Casassa91, Sukhorukov13, Pritchard12, Lishman13}.
157 gezelter 4243
158     Using molecular dynamics simulations, Samadashvili has recently shown
159     that when two smooth ice slabs slide past one another, a stable
160 gezelter 4254 liquid-like layer develops between them~\cite{Samadashvili13}. In a
161 gezelter 4245 previous study, our RNEMD simulations of ice-I$_\mathrm{h}$ shearing
162     through liquid water have provided quantitative estimates of the
163 gezelter 4254 solid-liquid kinetic friction coefficients~\cite{Louden13}. These
164 gezelter 4245 displayed a factor of two difference between the basal and prismatic
165     facets. The friction was found to be independent of shear direction
166     relative to the surface orientation. We attributed facet-based
167     difference in liquid-solid friction to the 6.5 \AA\ corrugation of the
168     prismatic face which reduces the effective surface area of the ice
169     that is in direct contact with liquid water.
170 gezelter 4243
171 gezelter 4245 In the sections that follow, we outline the methodology used to
172     simulate droplet-spreading dynamics using standard MD and tribological
173     properties using RNEMD simulations. These simulation methods give
174     complementary results that point to the prismatic and secondary prism
175     facets having roughly half of their surface area in direct contact
176     with the liquid.
177 gezelter 4243
178 gezelter 4245 \section{Methodology}
179     \subsection{Construction of the Ice / Water Interfaces}
180     To construct the four interfacial ice/water systems, a proton-ordered,
181     zero-dipole crystal of ice-I$_\mathrm{h}$ with exposed strips of
182 gezelter 4247 H-atoms was created using Structure 6 of Hirsch and Ojam\"{a}e's set
183     of orthorhombic representations for ice-I$_{h}$~\cite{Hirsch04}. This
184     crystal structure was cleaved along the four different facets. The
185     exposed face was reoriented normal to the $z$-axis of the simulation
186     cell, and the structures were and extended to form large exposed
187     facets in rectangular box geometries. Liquid water boxes were created
188     with identical dimensions (in $x$ and $y$) as the ice, with a $z$
189     dimension of three times that of the ice block, and a density
190     corresponding to 1 g / cm$^3$. Each of the ice slabs and water boxes
191     were independently equilibrated at a pressure of 1 atm, and the
192     resulting systems were merged by carving out any liquid water
193     molecules within 3 \AA\ of any atoms in the ice slabs. Each of the
194     combined ice/water systems were then equilibrated at 225K, which is
195     the liquid-ice coexistence temperature for SPC/E
196 gezelter 4254 water~\cite{Bryk02}. Reference \citealp{Louden13} contains a more
197     detailed explanation of the construction of similar ice/water
198     interfaces. The resulting dimensions as well as the number of ice and
199     liquid water molecules contained in each of these systems are shown in
200     Table \ref{tab:method}.
201 gezelter 4243
202 gezelter 4247 The SPC/E water model~\cite{Berendsen87} has been extensively
203 plouden 4246 characterized over a wide range of liquid
204 gezelter 4254 conditions~\cite{Arbuckle02,Kuang12}, and its phase diagram has been
205     well studied~\cite{Baez95,Bryk04b,Sanz04b,Fennell:2005fk}, With longer
206 gezelter 4247 cutoff radii and careful treatment of electrostatics, SPC/E mostly
207     avoids metastable crystalline morphologies like
208     ice-\textit{i}~\cite{Fennell:2005fk} and ice-B~\cite{Baez95}. The
209 gezelter 4254 free energies and melting
210     points~\cite{Baez95,Arbuckle02,Gay02,Bryk02,Bryk04b,Sanz04b,Fennell:2005fk,Fernandez06,Abascal07,Vrbka07}
211 gezelter 4247 of various other crystalline polymorphs have also been calculated.
212     Haymet \textit{et al.} have studied quiescent Ice-I$_\mathrm{h}$/water
213     interfaces using the SPC/E water model, and have seen structural and
214     dynamic measurements of the interfacial width that agree well with
215     more expensive water models, although the coexistence temperature for
216     SPC/E is still well below the experimental melting point of real
217     water~\cite{Bryk02}. Given the extensive data and speed of this model,
218     it is a reasonable choice even though the temperatures required are
219     somewhat lower than real ice / water interfaces.
220 gezelter 4245
221 gezelter 4247 \subsection{Droplet Simulations}
222     Ice interfaces with a thickness of $\sim$~20~\AA\ were created as
223 gezelter 4245 described above, but were not solvated in a liquid box. The crystals
224     were then replicated along the $x$ and $y$ axes (parallel to the
225 gezelter 4247 surface) until a large surface ($>$ 126 nm\textsuperscript{2}) had
226     been created. The sizes and numbers of molecules in each of the
227     surfaces is given in Table \ref{tab:method}. Weak translational
228 gezelter 4254 restraining potentials with spring constants of 1.5~$\mathrm{kcal\
229     mol}^{-1}\mathrm{~\AA}^{-2}$ (prismatic and pyramidal facets) or
230     4.0~$\mathrm{kcal\ mol}^{-1}\mathrm{~\AA}^{-2}$ (basal facet) were
231     applied to the centers of mass of each molecule in order to prevent
232     surface melting, although the molecules were allowed to reorient
233     freely. A water doplet containing 2048 SPC/E molecules was created
234     separately. Droplets of this size can produce agreement with the Young
235     contact angle extrapolated to an infinite drop size~\cite{Daub10}. The
236     surfaces and droplet were independently equilibrated to 225 K, at
237     which time the droplet was placed 3-5~\AA\ above the surface. Five
238     statistically independent simulations were carried out for each facet,
239     and the droplet was placed at unique $x$ and $y$ locations for each of
240     these simulations. Each simulation was 5~ns in length and was
241     conducted in the microcanonical (NVE) ensemble. Representative
242     configurations for the droplet on the prismatic facet are shown in
243     figure \ref{fig:Droplet}.
244 gezelter 4243
245 gezelter 4247 \subsection{Shearing Simulations (Interfaces in Bulk Water)}
246    
247     To perform the shearing simulations, the velocity shearing and scaling
248     variant of reverse non-equilibrium molecular dynamics (VSS-RNEMD) was
249     employed \cite{Kuang12}. This method performs a series of simultaneous
250     non-equilibrium exchanges of linear momentum and kinetic energy
251     between two physically-separated regions of the simulation cell. The
252     system responds to this unphysical flux with velocity and temperature
253     gradients. When VSS-RNEMD is applied to bulk liquids, transport
254     properties like the thermal conductivity and the shear viscosity are
255     easily extracted assuming a linear response between the flux and the
256     gradient. At the interfaces between dissimilar materials, the same
257     method can be used to extract \textit{interfacial} transport
258     properties (e.g. the interfacial thermal conductance and the
259     hydrodynamic slip length).
260    
261     The kinetic energy flux (producing a thermal gradient) is necessary
262     when performing shearing simulations at the ice-water interface in
263     order to prevent the frictional heating due to the shear from melting
264 gezelter 4254 the crystal. Reference \citealp{Louden13} provides more details on the
265     VSS-RNEMD method as applied to ice-water interfaces. A representative
266     configuration of the solvated prismatic facet being sheared through
267     liquid water is shown in figure \ref{fig:Shearing}.
268 gezelter 4247
269 gezelter 4254 The exchanges between the two regions were carried out every 2 fs
270     (e.g. every time step). This was done to minimize the magnitude of
271     each individual momentum exchange. Because individual VSS-RNEMD
272     exchanges conserve both total energy and linear momentum, the method
273     can be ``bolted-on'' to simulations in any ensemble. The simulations
274     of the pyramidal interface were performed under the canonical (NVT)
275     ensemble. When time correlation functions were computed, the RNEMD
276     simulations were done in the microcanonical (NVE) ensemble. All
277     simulations of the other interfaces were carried out in the
278     microcanonical ensemble.
279 gezelter 4247
280     \section{Results}
281     \subsection{Ice - Water Contact Angles}
282 plouden 4246
283     To determine the extent of wetting for each of the four crystal
284 gezelter 4247 facets, contact angles for liquid droplets on the ice surfaces were
285     computed using two methods. In the first method, the droplet is
286     assumed to form a spherical cap, and the contact angle is estimated
287     from the $z$-axis location of the droplet's center of mass
288     ($z_\mathrm{cm}$). This procedure was first described by Hautman and
289     Klein~\cite{Hautman91}, and was utilized by Hirvi and Pakkanen in
290     their investigation of water droplets on polyethylene and poly(vinyl
291     chloride) surfaces~\cite{Hirvi06}. For each stored configuration, the
292     contact angle, $\theta$, was found by inverting the expression for the
293     location of the droplet center of mass,
294 plouden 4246 \begin{equation}\label{contact_1}
295 gezelter 4247 \langle z_\mathrm{cm}\rangle = 2^{-4/3}R_{0}\bigg(\frac{1-cos\theta}{2+cos\theta}\bigg)^{1/3}\frac{3+cos\theta}{2+cos\theta} ,
296 plouden 4246 \end{equation}
297 gezelter 4247 where $R_{0}$ is the radius of the free water droplet.
298 plouden 4246
299 gezelter 4247 The second method for obtaining the contact angle was described by
300     Ruijter, Blake, and Coninck~\cite{Ruijter99}. This method uses a
301     cylindrical averaging of the droplet's density profile. A threshold
302     density of 0.5 g cm\textsuperscript{-3} is used to estimate the
303     location of the edge of the droplet. The $r$ and $z$-dependence of
304     the droplet's edge is then fit to a circle, and the contact angle is
305     computed from the intersection of the fit circle with the $z$-axis
306     location of the solid surface. Again, for each stored configuration,
307     the density profile in a set of annular shells was computed. Due to
308     large density fluctuations close to the ice, all shells located within
309     2 \AA\ of the ice surface were left out of the circular fits. The
310     height of the solid surface ($z_\mathrm{suface}$) along with the best
311 gezelter 4254 fitting origin ($z_\mathrm{droplet}$) and radius
312 gezelter 4247 ($r_\mathrm{droplet}$) of the droplet can then be used to compute the
313     contact angle,
314     \begin{equation}
315 gezelter 4254 \theta = 90 + \frac{180}{\pi} \sin^{-1}\left(\frac{z_\mathrm{droplet} -
316 gezelter 4247 z_\mathrm{surface}}{r_\mathrm{droplet}} \right).
317     \end{equation}
318     Both methods provided similar estimates of the dynamic contact angle,
319     although the first method is significantly less prone to noise, and
320     is the method used to report contact angles below.
321    
322     Because the initial droplet was placed above the surface, the initial
323 gezelter 4250 value of 180$^{\circ}$ decayed over time (See figure
324     \ref{fig:ContactAngle}). Each of these profiles were fit to a
325     biexponential decay, with a short-time contribution ($\tau_c$) that
326     describes the initial contact with the surface, a long time
327     contribution ($\tau_s$) that describes the spread of the droplet over
328     the surface, and a constant ($\theta_\infty$) to capture the
329     infinite-time estimate of the equilibrium contact angle,
330 gezelter 4247 \begin{equation}
331 gezelter 4250 \theta(t) = \theta_\infty + (180-\theta_\infty) \left[ a e^{-t/\tau_c} +
332     (1-a) e^{-t/\tau_s} \right]
333 gezelter 4247 \end{equation}
334 gezelter 4250 We have found that the rate for water droplet spreading across all
335     four crystal facets, $k_\mathrm{spread} = 1/\tau_s \approx$ 0.7
336     ns$^{-1}$. However, the basal and pyramidal facets produced estimated
337     equilibrium contact angles, $\theta_\infty \approx$ 35$^{o}$, while
338 gezelter 4247 prismatic and secondary prismatic had values for $\theta_\infty$ near
339 gezelter 4250 43$^{o}$ as seen in Table \ref{tab:kappa}.
340 gezelter 4247
341 gezelter 4254 These results indicate that the basal and pyramidal facets are more
342     hydrophilic than the prismatic and secondary prism facets, and
343     surprisingly, that the differential hydrophilicities of the crystal
344     facets is not reflected in the spreading rate of the droplet.
345 gezelter 4247
346 plouden 4246 % This is in good agreement with our calculations of friction
347     % coefficients, in which the basal
348     % and pyramidal had a higher coefficient of kinetic friction than the
349     % prismatic and secondary prismatic. Due to this, we beleive that the
350     % differences in friction coefficients can be attributed to the varying
351     % hydrophilicities of the facets.
352    
353     \subsection{Coefficient of friction of the interfaces}
354 gezelter 4254 In a bulk fluid, the shear viscosity, $\eta$, can be determined
355     assuming a linear response to a shear stress,
356 plouden 4246 \begin{equation}\label{Shenyu-11}
357 gezelter 4254 j_{z}(p_{x}) = \eta \frac{\partial v_{x}}{\partial z}.
358 plouden 4246 \end{equation}
359 gezelter 4254 Here $j_{z}(p_{x})$ is the flux (in $x$-momentum) that is transferred
360     in the $z$ direction (i.e. the shear stress). The RNEMD simulations
361     impose an artificial momentum flux between two regions of the
362     simulation, and the velocity gradient is the fluid's response. This
363     technique has now been applied quite widely to determine the
364     viscosities of a number of bulk fluids~\cite{}.
365    
366     At the interface between two phases (e.g. liquid / solid) the same
367     momentum flux creates a velocity difference between the two materials,
368     and this can be used to define an interfacial friction coefficient
369     ($\kappa$),
370     \begin{equation}\label{Shenyu-13}
371     j_{z}(p_{x}) = \kappa \left[ v_{x}(liquid) - v_{x}(solid) \right]
372 plouden 4246 \end{equation}
373 gezelter 4254 where $v_{x}(liquid)$ and $v_{x}(solid)$ are the velocities measured
374     directly adjacent to the interface.
375    
376     The simulations described here contain significant quantities of both
377     liquid and solid phases, and the momentum flux must traverse a region
378     of the liquid that is simultaneously under a thermal gradient. Since
379     the liquid has a temperature-dependent shear viscosity, $\eta(T)$,
380     estimates of the solid-liquid friction coefficient can be obtained if
381     one knows the viscosity of the liquid at the interface (i.e. at the
382     melting temperature, $T_m$),
383 plouden 4246 \begin{equation}\label{kappa-2}
384 gezelter 4254 \kappa = \frac{\eta(T_{m})}{\left[v_{x}(fluid)-v_{x}(solid)\right]}\left(\frac{\partial v_{x}}{\partial z}\right).
385 plouden 4246 \end{equation}
386 gezelter 4254 For SPC/E, the melting temperature of Ice-I$_\mathrm{h}$ is estimated
387     to be 225~K~\cite{Bryk02}. To obtain the value of
388     $\eta(225\mathrm{~K})$ for the SPC/E model, a $31.09 \times 29.38
389     \times 124.39$ \AA\ box with 3744 water molecules in a disordered
390     configuration was equilibrated to 225~K, and five
391     statistically-independent shearing simulations were performed (with
392 plouden 4255 imposed fluxes that spanned a range of $3 \rightarrow 13 \mathrm{~m~s}^{-1}$ ). Each simulation
393 gezelter 4254 was conducted in the microcanonical ensemble with total simulation
394     times of 5 ns. The VSS-RNEMD exchanges were carried out every 2 fs. We
395     estimate $\eta(225\mathrm{~K})$ to be 0.0148 $\pm$ 0.0007 Pa s for
396     SPC/E, roughly ten times larger than the shear viscosity previously
397     computed at 280~K~\cite{Kuang12}.
398 plouden 4246
399 gezelter 4254 The interfacial friction coefficient, $\kappa$, can equivalently be
400     expressed as the ratio of the viscosity of the fluid to the
401     hydrodynamic slip length, $\delta$, which is an indication of strength
402     of the interactions between the solid and liquid phases,
403 plouden 4246 \begin{equation}\label{kappa-3}
404     \kappa = \frac{\eta}{\delta}
405     \end{equation}
406 gezelter 4254 The connection between slip length and surface hydrophobicity is not
407     yet clear. In some simulations, the slip length has been found to have
408     a link to the effective surface hydrophobicity~\cite{Sendner:2009uq},
409     although Ho \textit{et al.} have found that liquid water can also slip
410     on hydrophilic surfaces~\cite{Ho:2011zr}. Experimental evidence for a
411     direct tie between slip length and hydrophobicity is also not
412     definitive. Total-internal reflection particle image velocimetry
413     (TIR-PIV) studies have suggested that there is a link between slip
414     length and effective
415     hydrophobicity~\cite{Lasne:2008vn,Bouzigues:2008ys}. However, recent
416     surface sensitive cross-correlation spectroscopy (TIR-FCCS)
417     measurements have seen similar slip behavior for both hydrophobic and
418     hydrophilic surfaces~\cite{Schaeffel:2013kx}.
419 plouden 4246
420 gezelter 4254 In each of the systems studied here, the interfacial temperature was
421     kept fixed to 225K, which ensured the viscosity of the fluid at the
422     interace was identical. Thus, any significant variation in $\kappa$
423     between the systems is a direct indicator of the slip length and the
424     effective interaction strength between the solid and liquid phases.
425 gezelter 4243
426 gezelter 4254 The calculated $\kappa$ values found for the four crystal facets of
427     Ice-I$_\mathrm{h}$ are shown in Table \ref{tab:kappa}. The basal and
428     pyramidal facets were found to have similar values of $\kappa \approx
429     6$ ($\times 10^{-4} \mathrm{amu~\AA}^{-2}\mathrm{fs}^{-1}$), while the
430     prismatic and secondary prism facets exhibited $\kappa \approx 3$
431     ($\times 10^{-4} \mathrm{amu~\AA}^{-2}\mathrm{fs}^{-1}$). These
432     results are also essentially independent of shearing direction
433     relative to features on the surface of the facets. The friction
434     coefficients indicate that the basal and pyramidal facets have
435     significantly stronger interactions with liquid water than either of
436     the two prismatic facets. This is in agreement with the contact angle
437     results above - both of the high-friction facets exhbited smaller
438     contact angles, suggesting that the solid-liquid friction is
439     correlated with the hydrophilicity of these facets.
440 gezelter 4243
441 gezelter 4254 \subsection{Structural measures of interfacial width under shear}
442     One of the open questions about ice/water interfaces is whether the
443     thickness of the 'slush-like' quasi-liquid layer (QLL) depends on the
444     facet of ice presented to the water. In the QLL region, the water
445     molecules are ordered differently than in either the solid or liquid
446     phases, and also exhibit distinct dynamical behavior. The width of
447     this quasi-liquid layer has been estimated by finding the distance
448     over which structural order parameters or dynamic properties change
449     from their bulk liquid values to those of the solid ice. The
450     properties used to find interfacial widths have included the local
451     density, the diffusion constant, and the translational and
452     orientational order
453     parameters~\cite{Karim88,Karim90,Hayward01,Hayward02,Bryk02,Gay02,Louden13}.
454    
455     The VSS-RNEMD simulations impose thermal and velocity gradients.
456     These gradients perturb the momenta of the water molecules, so
457     parameters that depend on translational motion are often measuring the
458     momentum exchange, and not physical properties of the interface. As a
459     structural measure of the interface, we have used the local
460     tetrahedral order parameter to estimate the width of the interface.
461     This quantity was originally described by Errington and
462     Debenedetti~\cite{Errington01} and has been used in bulk simulations
463     by Kumar \textit{et al.}~\cite{Kumar09}. It has previously been used
464     in ice/water interfaces by by Bryk and Haymet~\cite{Bryk04b}.
465    
466     To determine the structural widths of the interfaces under shear, each
467     of the systems was divided into 100 bins along the $z$-dimension, and
468     the local tetrahedral order parameter (Eq. 5 in Reference
469     \citealp{Louden13}) was time-averaged in each bin for the duration of
470     the shearing simulation. The spatial dependence of this order
471     parameter, $q(z)$, is the tetrahedrality profile of the interface. A
472     representative profile for the pyramidal facet is shown in circles in
473     panel $a$ of figure \ref{fig:pyrComic}. The $q(z)$ function has a
474     range of $(0,1)$, where a value of unity indicates a perfectly
475     tetrahedral environment. The $q(z)$ for the bulk liquid was found to
476     be $\approx~0.77$, while values of $\approx~0.92$ were more common in
477     the ice. The tetrahedrality profiles were fit using a hyperbolic
478     tangent function (see Eq. 6 in Reference \citealp{Louden13}) designed
479     to smoothly fit the bulk to ice transition while accounting for the
480     weak thermal gradient. In panels $b$ and $c$, the resulting thermal
481     and velocity gradients from an imposed kinetic energy and momentum
482     fluxes can be seen. The vertical dotted lines traversing all three
483     panels indicate the midpoints of the interface as determined by the
484     tetrahedrality profiles.
485 gezelter 4243
486 gezelter 4254 We find the interfacial width to be $3.2 \pm 0.2$ \AA\ (pyramidal) and
487     $3.2 \pm 0.2$ \AA\ (secondary prism) for the control systems with no
488     applied momentum flux. This is similar to our previous results for the
489     interfacial widths of the quiescent basal ($3.2 \pm 0.4$ \AA) and
490     prismatic systems ($3.6 \pm 0.2$ \AA).
491 gezelter 4243
492 gezelter 4254 Over the range of shear rates investigated, $0.4 \rightarrow
493     6.0~\mathrm{~m~s}^{-1}$ for the pyramidal system and $0.6 \rightarrow
494     5.2~\mathrm{~m~s}^{-1}$ for the secondary prism, we found no
495     significant change in the interfacial width. The mean interfacial
496     widths are collected in table \ref{tab:kappa}. This follows our
497     previous findings of the basal and prismatic systems, in which the
498     interfacial widths of the basal and prismatic facets were also found
499     to be insensitive to the shear rate~\cite{Louden13}.
500 gezelter 4243
501 gezelter 4254 The similarity of these interfacial width estimates indicate that the
502     particular facet of the exposed ice crystal has little to no effect on
503     how far into the bulk the ice-like structural ordering persists. Also,
504     it appears that for the shearing rates imposed in this study, the
505     interfacial width is not structurally modified by the movement of
506     water over the ice.
507 gezelter 4243
508 gezelter 4254 \subsection{Dynamic measures of interfacial width under shear}
509 gezelter 4243 The orientational time correlation function,
510     \begin{equation}\label{C(t)1}
511     C_{2}(t)=\langle P_{2}(\mathbf{u}(0)\cdot \mathbf{u}(t))\rangle,
512     \end{equation}
513     helps indicate the local environment around the water molecules. The function
514     begins with an initial value of unity, and decays to zero as the water molecule
515     loses memory of its former orientation. Observing the rate at which this decay
516     occurs can provide insight to the mechanism and timescales for the relaxation.
517     In eq. \eqref{C(t)1}, $P_{2}$ is the second-order Legendre polynomial, and
518     $\mathbf{u}$ is the bisecting HOH vector. The angle brackets indicate
519     an ensemble average over all the water molecules in a given spatial region.
520    
521     To investigate the dynamics of the water molecules across the interface, the
522     systems were divided in the $z$-dimension into bins, each $\approx$ 3 \AA\
523     wide, and eq. \eqref{C(t)1} was computed for each of the bins. A water
524     molecule was allocated to a particular bin if it was initially in the bin
525     at time zero. To compute eq. \eqref{C(t)1}, each 0.5 ns simulation was
526     followed by an additional 200 ps NVE simulation during which the
527     position and orientations of each molecule were recorded every 0.1 ps.
528    
529     The data obtained for each bin was then fit to a triexponential decay
530     with the three decay constants
531     $\tau_{short}$ corresponding to the librational motion of the water
532     molecules, $\tau_{middle}$ corresponding to jumps between the breaking and
533     making of hydrogen bonds, and $\tau_{long}$ corresponding to the translational
534     motion of the water molecules. An additive constant in the fit accounts
535     for the water molecules trapped in the ice which do not experience any
536     long-time orientational decay.
537    
538 gezelter 4250 In Figure \ref{fig:PyrOrient} we see the $z$-coordinate profiles for
539     the three decay constants, $\tau_{short}$ (panel a), $\tau_{middle}$
540     (panel b), and $\tau_{long}$ (panel c) for the pyramidal and secondary
541     prismatic systems respectively. The control experiments (no shear) are
542     shown with circles, and an experiment with an imposed momentum flux is
543     shown with squares. The vertical dotted line traversing all three
544     panels denotes the midpoint of the interface determined using the
545     local tetrahedral order parameter. In the liquid regions of both
546     systems, we see that $\tau_{middle}$ and $\tau_{long}$ have
547     approximately consistent values of $3-6$ ps and $30-40$ ps,
548     resepctively, and increase in value as we approach the
549     interface. Conversely, in panel a, we see that $\tau_{short}$
550     decreases from the liquid value of $72-76$ fs as we approach the
551     interface. We believe this speed up is due to the constrained motion
552     of librations closer to the interface. Both the approximate values for
553     the decays and trends approaching the interface match those reported
554     previously for the basal and prismatic interfaces.
555 gezelter 4243
556     As done previously, we have attempted to quantify the distance, $d_{pyramidal}$
557     and $d_{secondary prismatic}$, from the
558     interface that the deviations from the bulk liquid values begin. This was done
559     by fitting the orientational decay constant $z$-profiles by
560     \begin{equation}\label{tauFit}
561     \tau(z)\approx\tau_{liquid}+(\tau_{wall}-\tau_{liquid})e^{-(z-z_{wall})/d}
562     \end{equation}
563     where $\tau_{liquid}$ and $\tau_{wall}$ are the liquid and projected wall
564     values of the decay constants, $z_{wall}$ is the location of the interface,
565     and $d$ is the displacement from the interface at which these deviations
566     occur. The values for $d_{pyramidal}$ and $d_{secondary prismatic}$ were
567     determined
568     for each of the decay constants, and then averaged for better statistics
569     ($\tau_{middle}$ was ommitted for secondary prismatic). For the pyramidal
570     system,
571     $d_{pyramidal}$ was found to be 2.7 \AA\ for both the control and the sheared
572     system. We found $d_{secondary prismatic}$ to be slightly larger than
573     $d_{pyramidal}$ for both the control and with an applied shear, with
574     displacements of $4$ \AA\ for the control system and $3$ \AA\ for the
575     experiment with the imposed momentum flux. These values are consistent with
576     those found for the basal ($d_{basal}\approx2.9$ \AA\ ) and prismatic
577     ($d_{prismatic}\approx3.5$ \AA\ ) systems.
578    
579    
580     \section{Conclusion}
581     We present the results of molecular dynamics simulations of the basal,
582     prismatic, pyrmaidal
583     and secondary prismatic facets of an SPC/E model of the
584     Ice-I$_\mathrm{h}$/water interface, and show that the differential
585     coefficients of friction among the four facets are due to their
586     relative hydrophilicities by means
587     of water contact angle calculations. To obtain the coefficients of
588     friction, the ice was sheared through the liquid
589     water while being exposed to a thermal gradient to maintain a stable
590     interface by using the minimally perturbing VSS RNEMD method. Water
591     contact angles are obtained by fitting the spreading of a liquid water
592     droplet over the crystal facets.
593    
594     In agreement with our previous findings for the basal and prismatic facets, the interfacial
595     width of the prismatic and secondary prismatic crystal faces were
596     found to be independent of shear rate as measured by the local
597     tetrahedral order parameter. This width was found to be
598     3.2~$\pm$ 0.2~\AA\ for both the pyramidal and the secondary prismatic systems.
599     These values are in good agreement with our previously calculated interfacial
600     widths for the basal (3.2~$\pm$ 0.4~\AA\ ) and prismatic (3.6~$\pm$ 0.2~\AA\ )
601     systems.
602    
603     Orientational dynamics of the Ice-I$_\mathrm{h}$/water interfaces were studied
604     by calculation of the orientational time correlation function at varying
605     displacements normal to the interface. The decays were fit
606     to a tri-exponential decay, where the three decay constants correspond to
607     the librational motion of the molecules driven by the restoring forces of
608     existing hydrogen bonds ($\tau_{short}$ $\mathcal{O}$(10 fs)), jumps between
609     two different hydrogen bonds ($\tau_{middle}$ $\mathcal{O}$(1 ps)), and
610     translational motion of the molecules ($\tau_{long}$ $\mathcal{O}$(100 ps)).
611     $\tau_{short}$ was found to decrease approaching the interface due to the
612     constrained motion of the molecules as the local environment becomes more
613     ice-like. Conversely, the two longer-time decay constants were found to
614     increase at small displacements from the interface. As seen in our previous
615     work on the basal and prismatic facets, there appears to be a dynamic
616     interface width at which deviations from the bulk liquid values occur.
617     We had previously found $d_{basal}$ and $d_{prismatic}$ to be approximately
618     2.8~\AA\ and 3.5~\AA. We found good agreement of these values for the
619     pyramidal and secondary prismatic systems with $d_{pyramidal}$ and
620     $d_{secondary prismatic}$ to be 2.7~\AA\ and 3~\AA. For all four of the
621     facets, no apparent dependence of the dynamic width on the shear rate was
622     found.
623    
624     The interfacial friction coefficient, $\kappa$, was determined for each facet
625     interface. We were able to reach an expression for $\kappa$ as a function of
626     the velocity profile of the system which is scaled by the viscosity of the liquid
627     at 225 K. In doing so, we have obtained an expression for $\kappa$ which is
628     independent of temperature differences of the liquid water at far displacements
629     from the interface. We found the basal and pyramidal facets to have
630     similar $\kappa$ values, of $\kappa \approx$ 6
631     (x$10^{-4}$amu \AA\textsuperscript{-2} fs\textsuperscript{-1}). However, the
632     prismatic and secondary prismatic facets were found to have $\kappa$ values of
633     $\kappa \approx$ 3 (x$10^{-4}$amu \AA\textsuperscript{-2} fs\textsuperscript{-1}).
634     Believing this difference was due to the relative hydrophilicities of
635     the crystal faces, we have calculated the infinite decay of the water
636     contact angle, $\theta_{\infty}$, by watching the spreading of a water
637     droplet over the surface of the crystal facets. We have found
638     $\theta_{\infty}$ for the basal and pyramidal faces to be $\approx$ 34
639     degrees, while obtaining $\theta_{\infty}$ of $\approx$ 43 degrees for
640     the prismatic and secondary prismatic faces. This indicates that the
641     basal and pyramidal faces of ice-I$_\mathrm{h}$ are more hydrophilic
642     than the prismatic and secondary prismatic. These results also seem to
643     explain the differential friction coefficients obtained through the
644     shearing simulations, namely, that the coefficients of friction of the
645     ice-I$_\mathrm{h}$ crystal facets are governed by their inherent
646     hydrophilicities.
647    
648    
649     \begin{acknowledgments}
650     Support for this project was provided by the National
651     Science Foundation under grant CHE-1362211. Computational time was
652     provided by the Center for Research Computing (CRC) at the
653     University of Notre Dame.
654     \end{acknowledgments}
655    
656     \bibliography{iceWater}
657     % *****************************************
658     % There is significant interest in the properties of ice/ice and ice/water
659     % interfaces in the geophysics community. Most commonly, the results of shearing
660     % two ice blocks past one
661     % another\cite{Casassa91, Sukhorukov13, Pritchard12, Lishman13} or the shearing
662     % of ice through water\cite{Cuffey99, Bell08}. Using molecular dynamics
663     % simulations, Samadashvili has recently shown that when two smooth ice slabs
664     % slide past one another, a stable liquid-like layer develops between
665     % them\cite{Samadashvili13}. To fundamentally understand these processes, a
666     % molecular understanding of the ice/water interfaces is needed.
667    
668     % Investigation of the ice/water interface is also crucial in understanding
669     % processes such as nucleation, crystal
670     % growth,\cite{Han92, Granasy95, Vanfleet95} and crystal
671     % melting\cite{Weber83, Han92, Sakai96, Sakai96B}. Insight gained to these
672     % properties can also be applied to biological systems of interest, such as
673     % the behavior of the antifreeze protein found in winter
674     % flounder,\cite{Wierzbicki07,Chapsky97} and certain terrestial
675     % arthropods.\cite{Duman:2001qy,Meister29012013} Elucidating the properties which
676     % give rise to these processes through experimental techniques can be expensive,
677     % complicated, and sometimes infeasible. However, through the use of molecular
678     % dynamics simulations much of the problems of investigating these properties
679     % are alleviated.
680    
681     % Understanding ice/water interfaces inherently begins with the isolated
682     % systems. There has been extensive work parameterizing models for liquid water,
683     % such as the SPC\cite{Berendsen81}, SPC/E\cite{Berendsen87},
684     % TIP4P\cite{Jorgensen85}, TIP4P/2005\cite{Abascal05},
685     % ($\dots$), and more recently, models for simulating
686     % the solid phases of water, such as the TIP4P/Ice\cite{Abascal05b} model. The
687     % melting point of various crystal structures of ice have been calculated for
688     % many of these models
689     % (SPC\cite{Karim90,Abascal07}, SPC/E\cite{Baez95,Arbuckle02,Gay02,Bryk02,Bryk04b,Sanz04b,Gernandez06,Abascal07,Vrbka07}, TIP4P\cite{Karim88,Gao00,Sanz04,Sanz04b,Koyama04,Wang05,Fernandez06,Abascal07}, TIP5P\cite{Sanz04,Koyama04,Wang05,Fernandez06,Abascal07}),
690     % and the partial or complete phase diagram for the model has been determined
691     % (SPC/E\cite{Baez95,Bryk04b,Sanz04b}, TIP4P\cite{Sanz04,Sanz04b,Koyama04}, TIP5P\cite{Sanz04,Koyama04}).
692     % Knowing the behavior and melting point for these models has enabled an initial
693     % investigation of ice/water interfaces.
694    
695     % The Ice-I$_\mathrm{h}$/water quiescent interface has been extensively studied
696     % over the past 30 years by theory and experiment. Haymet \emph{et al.} have
697     % done significant work characterizing and quantifying the width of these
698     % interfaces for the SPC,\cite{Karim90} SPC/E,\cite{Gay02,Bryk02},
699     % CF1,\cite{Hayward01,Hayward02} and TIP4P\cite{Karim88} models for water. In
700     % recent years, Haymet has focused on investigating the effects cations and
701     % anions have on crystal nucleaion and
702     % melting.\cite{Bryk04,Smith05,Wilson08,Wilson10} Nada and Furukawa have studied
703     % the the basal- and prismatic-water interface width\cite{Nada95}, crystal
704     % surface restructuring at temperatures approaching the melting
705     % point\cite{Nada00}, and the mechanism of ice growth inhibition by antifreeze
706     % proteins\cite{Nada08,Nada11,Nada12}. Nada has developed a six-site water model
707     % for ice/water interfaces near the melting point\cite{Nada03}, and studied the
708     % dependence of crystal growth shape on applied pressure\cite{Nada11b}. Using
709     % this model, Nada and Furukawa have established differential
710     % growth rates for the basal, prismatic, and secondary prismatic facets of
711     % Ice-I$_\mathrm{h}$ and found their origins due to a reordering of the hydrogen
712     % bond network in water near the interface\cite{Nada05}. While the work
713     % described so far has mainly focused on bulk water on ice, there is significant
714     % interest in thin films of water on ice surfaces as well.
715    
716     % It is well known that the surface of ice exhibits a premelting layer at
717     % temperatures near the melting point, often called a quasi-liquid layer (QLL).
718     % Molecular dynamics simulations of the facets of ice-I$_\mathrm{h}$ exposed
719     % to vacuum performed by Conde, Vega and Patrykiejew have found QLL widths of
720     % approximately 10 \AA\ at 3 K below the melting point\cite{Conde08}.
721     % Similarly, Limmer and Chandler have used course grain simulations and
722     % statistical field theory to estimated QLL widths at the same temperature to
723     % be about 3 nm\cite{Limmer14}.
724     % Recently, Sazaki and Furukawa have developed an experimental technique with
725     % sufficient spatial and temporal resolution to visulaize and quantitatively
726     % analyze QLLs on ice crystals at temperatures near melting\cite{Sazaki10}. They
727     % have found the width of the QLLs perpindicular to the surface at -2.2$^{o}$C
728     % to be $\mathcal{O}$(\AA). They have also seen the formation of two immiscible
729     % QLLs, which displayed different stabilities and dynamics on the crystal
730     % surface\cite{Sazaki12}. Knowledge of the hydrophilicities of each
731     % of the crystal facets would help further our understanding of the properties
732     % and dynamics of the QLLs.
733    
734     % Presented here is the follow up to our previous paper\cite{Louden13}, in which
735     % the basal and prismatic facets of an ice-I$_\mathrm{h}$/water interface were
736     % investigated where the ice was sheared relative to the liquid. By using a
737     % recently developed velocity shearing and scaling approach to reverse
738     % non-equilibrium molecular dynamics (VSS-RNEMD), simultaneous temperature and
739     % velocity gradients can be applied to the system, which allows for measurment
740     % of friction and thermal transport properties while maintaining a stable
741     % interfacial temperature\cite{Kuang12}. Structural analysis and dynamic
742     % correlation functions were used to probe the interfacial response to a shear,
743     % and the resulting solid/liquid kinetic friction coefficients were reported.
744     % In this paper we present the same analysis for the pyramidal and secondary
745     % prismatic facets, and show that the differential interfacial friction
746     % coefficients for the four facets of ice-I$_\mathrm{h}$ are determined by their
747     % relative hydrophilicity by means of dynamics water contact angle
748     % simulations.
749    
750     % The local tetrahedral order parameter, $q(z)$, is given by
751     % \begin{equation}
752     % q(z) = \int_0^L \sum_{k=1}^{N} \Bigg(1 -\frac{3}{8}\sum_{i=1}^{3}
753     % \sum_{j=i+1}^{4} \bigg(\cos\psi_{ikj}+\frac{1}{3}\bigg)^2\Bigg)
754     % \delta(z_{k}-z)\mathrm{d}z \Bigg/ N_z
755     % \label{eq:qz}
756     % \end{equation}
757     % where $\psi_{ikj}$ is the angle formed between the oxygen sites of molecules
758     % $i$,$k$, and $j$, where the centeral oxygen is located within molecule $k$ and
759     % molecules $i$ and $j$ are two of the closest four water molecules
760     % around molecule $k$. All four closest neighbors of molecule $k$ are also
761     % required to reside within the first peak of the pair distribution function
762     % for molecule $k$ (typically $<$ 3.41 \AA\ for water).
763     % $N_z = \int\delta(z_k - z) \mathrm{d}z$ is a normalization factor to account
764     % for the varying population of molecules within each finite-width bin.
765    
766    
767     % The hydrophobicity or hydrophilicity of a surface can be described by the
768     % extent a droplet of water wets the surface. The contact angle formed between
769     % the solid and the liquid, $\theta$, which relates the free energies of the
770     % three interfaces involved, is given by Young's equation.
771     % \begin{equation}\label{young}
772     % \cos\theta = (\gamma_{sv} - \gamma_{sl})/\gamma_{lv}
773     % \end{equation}
774     % Here $\gamma_{sv}$, $\gamma_{sl}$, and $\gamma_{lv}$ are the free energies
775     % of the solid/vapor, solid/liquid, and liquid/vapor interfaces respectively.
776     % Large contact angles ($\theta$ $\gg$ 90\textsuperscript{o}) correspond to low
777     % wettability and hydrophobic surfaces, while small contact angles
778     % ($\theta$ $\ll$ 90\textsuperscript{o}) correspond to high wettability and
779     % hydrophilic surfaces. Experimentally, measurements of the contact angle
780     % of sessile drops has been used to quantify the extent of wetting on surfaces
781     % with thermally selective wetting charactaristics\cite{Tadanaga00,Liu04,Sun04},
782     % as well as nano-pillared surfaces with electrically tunable Cassie-Baxter and
783     % Wenzel states\cite{Herbertson06,Dhindsa06,Verplanck07,Ahuja08,Manukyan11}.
784     % Luzar and coworkers have done significant work modeling these transitions on
785     % nano-patterned surfaces\cite{Daub07,Daub10,Daub11,Ritchie12}, and have found
786     % the change in contact angle to be due to the external field perturbing the
787     % hydrogen bonding of the liquid/vapor interface\cite{Daub07}.
788    
789    
790    
791     \end{article}
792    
793     \begin{figure}
794 gezelter 4247 \includegraphics[width=\linewidth]{Droplet}
795     \caption{\label{fig:Droplet} Computational model of a droplet of
796     liquid water spreading over the prismatic $\{1~0~\bar{1}~0\}$ facet
797 gezelter 4250 of ice, before (left) and 2.6 ns after (right) being introduced to the
798 gezelter 4247 surface. The contact angle ($\theta$) shrinks as the simulation
799     proceeds, and the long-time behavior of this angle is used to
800     estimate the hydrophilicity of the facet.}
801     \end{figure}
802    
803     \begin{figure}
804 gezelter 4250 \includegraphics[width=2in]{Shearing}
805 gezelter 4247 \caption{\label{fig:Shearing} Computational model of a slab of ice
806 gezelter 4250 being sheared through liquid water. In this figure, the ice is
807     presenting two copies of the prismatic $\{1~0~\bar{1}~0\}$ facet
808     towards the liquid phase. The RNEMD simulation exchanges both
809     linear momentum (indicated with arrows) and kinetic energy between
810     the central box and the box that spans the cell boundary. The
811     system responds with weak thermal gradient and a velocity profile
812     that shears the ice relative to the surrounding liquid.}
813 gezelter 4247 \end{figure}
814    
815     \begin{figure}
816 gezelter 4250 \includegraphics[width=\linewidth]{ContactAngle}
817     \caption{\label{fig:ContactAngle} The dynamic contact angle of a
818     droplet after approaching each of the four ice facets. The decay to
819     an equilibrium contact angle displays similar dynamics. Although
820     all the surfaces are hydrophilic, the long-time behavior stabilizes
821     to significantly flatter droplets for the basal and pyramidal
822     facets. This suggests a difference in hydrophilicity for these
823     facets compared with the two prismatic facets.}
824 gezelter 4243 \end{figure}
825    
826     \begin{figure}
827 gezelter 4250 \includegraphics[width=\linewidth]{Pyr_comic_strip}
828     \caption{\label{fig:pyrComic} Properties of the pyramidal interface
829     being sheared through water at 3.8 ms\textsuperscript{-1}. Lower
830     panel: the local tetrahedral order parameter, $q(z)$, (circles) and
831     the hyperbolic tangent fit (turquoise line). Middle panel: the
832     imposed thermal gradient required to maintain a fixed interfacial
833     temperature of 225 K. Upper panel: the transverse velocity gradient
834     that develops in response to an imposed momentum flux. The vertical
835     dotted lines indicate the locations of the midpoints of the two
836     interfaces.}
837 gezelter 4243 \end{figure}
838    
839 gezelter 4250 % \begin{figure}
840     % \includegraphics[width=\linewidth]{SP_comic_strip}
841     % \caption{\label{fig:spComic} The secondary prismatic interface with a shear
842     % rate of 3.5 \
843     % ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:pyrComic}.}
844     % \end{figure}
845    
846 gezelter 4243 \begin{figure}
847     \includegraphics[width=\linewidth]{Pyr-orient}
848 gezelter 4250 \caption{\label{fig:PyrOrient} The three decay constants of the
849     orientational time correlation function, $C_2(t)$, for water as a
850     function of distance from the center of the ice slab. The vertical
851     dashed line indicates the edge of the pyramidal ice slab determined
852     by the local order tetrahedral parameter. The control (circles) and
853     sheared (squares) simulations were fit using shifted-exponential
854     decay (see Eq. 9 in Ref. \citealp{Louden13}).}
855 gezelter 4243 \end{figure}
856    
857 gezelter 4250 % \begin{figure}
858     % \includegraphics[width=\linewidth]{SP-orient-less}
859     % \caption{\label{fig:SPorient} Decay constants for $C_2(t)$ at the secondary
860     % prismatic face. Panel descriptions match those in \ref{fig:PyrOrient}.}
861     % \end{figure}
862 gezelter 4243
863    
864     \begin{table}[h]
865     \centering
866 gezelter 4247 \caption{Sizes of the droplet and shearing simulations. Cell
867     dimensions are measured in \AA. \label{tab:method}}
868     \begin{tabular}{r|cccc|ccccc}
869     \toprule
870     \multirow{2}{*}{Interface} & \multicolumn{4}{c|}{Droplet} & \multicolumn{5}{c}{Shearing} \\
871     & $N_\mathrm{ice}$ & $N_\mathrm{droplet}$ & $L_x$ & $L_y$ & $N_\mathrm{ice}$ & $N_\mathrm{liquid}$ & $L_x$ & $L_y$ & $L_z$ \\
872     \midrule
873     Basal $\{0001\}$ & 12960 & 2048 & 134.70 & 140.04 & 900 & 1846 & 23.87 & 35.83 & 98.64 \\
874     Pyramidal $\{2~0~\bar{2}~1\}$ & 11136 & 2048 & 143.75 & 121.41 & 1216 & 2203 & 37.47 & 29.50 & 93.02 \\
875     Prismatic $\{1~0~\bar{1}~0\}$ & 9900 & 2048 & 110.04 & 115.00 & 3000 & 5464 & 35.95 & 35.65 & 205.77 \\
876     Secondary Prism $\{1~1~\bar{2}~0\}$ & 11520 & 2048 & 146.72 & 124.48 & 3840 & 8176 & 71.87 & 31.66 & 161.55 \\
877     \bottomrule
878 gezelter 4243 \end{tabular}
879     \end{table}
880    
881    
882     \begin{table}[h]
883     \centering
884 gezelter 4247 \caption{Structural and dynamic properties of the interfaces of Ice-I$_\mathrm{h}$
885     with water. Liquid-solid friction coefficients ($\kappa_x$ and $\kappa_y$) are expressed in 10\textsuperscript{-4} amu
886     \AA\textsuperscript{-2} fs\textsuperscript{-1}. \label{tab:kappa}}
887     \begin{tabular}{r|cc|cccc}
888     \toprule
889     \multirow{2}{*}{Interface} & \multicolumn{2}{c|}{Droplet} & \multicolumn{4}{c}{Shearing}\\
890     & $\theta_{\infty}$ ($^\circ$) & $k_\mathrm{spread}$ (ns\textsuperscript{-1}) &
891     $\kappa_{x}$ & $\kappa_{y}$ & $d_{q_{z}}$ (\AA) & $d_{\tau}$ (\AA) \\
892     \midrule
893 plouden 4251 Basal $\{0001\}$ & $34.1 \pm 0.9$ &$0.60 \pm 0.07$
894 plouden 4253 & $5.9 \pm 0.3$ & $6.5 \pm 0.8$ & $3.2 \pm 0.4$ & $2 \pm 1$ \\
895 gezelter 4250 Pyramidal $\{2~0~\bar{2}~1\}$ & $35 \pm 3$ & $0.7 \pm 0.1$ &
896 plouden 4253 $5.8 \pm 0.4$ & $6.1 \pm 0.5$ & $3.2 \pm 0.2$ & $2.5 \pm 0.3$\\
897 plouden 4251 Prismatic $\{1~0~\bar{1}~0\}$ & $45 \pm 3$ & $0.75 \pm 0.09$ &
898     $3.0 \pm 0.2$ & $3.0 \pm 0.1$ & $3.6 \pm 0.2$ & $4 \pm 2$ \\
899 plouden 4252 Secondary Prism $\{1~1~\bar{2}~0\}$ & $43 \pm 2$ & $0.69 \pm 0.03$ &
900 plouden 4253 $3.5 \pm 0.1$ & $3.3 \pm 0.2$ & $3.2 \pm 0.2$ & $5 \pm 3$ \\
901 gezelter 4247 \bottomrule
902 gezelter 4243 \end{tabular}
903     \end{table}
904    
905     \end{document}