ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceWater/iceWater.tex
Revision: 3970
Committed: Thu Oct 24 15:00:57 2013 UTC (11 years, 9 months ago) by gezelter
Content type: application/x-tex
File size: 41054 byte(s)
Log Message:
Latest version

File Contents

# Content
1 \documentclass[11pt]{article}
2 \usepackage{amsmath}
3 \usepackage{amssymb}
4 \usepackage{setspace}
5 %\usepackage{endfloat}
6 \usepackage{caption}
7 %\usepackage{epsf}
8 %\usepackage{tabularx}
9 \usepackage{graphicx}
10 \usepackage{multirow}
11 \usepackage{wrapfig}
12 %\usepackage{booktabs}
13 %\usepackage{bibentry}
14 %\usepackage{mathrsfs}
15 %\usepackage[ref]{overcite}
16 \usepackage[square, comma, sort&compress]{natbib}
17 \usepackage{url}
18 \pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
19 \evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
20 9.0in \textwidth 6.5in \brokenpenalty=10000
21
22 % double space list of tables and figures
23 %\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}}
24 \setlength{\abovecaptionskip}{20 pt}
25 \setlength{\belowcaptionskip}{30 pt}
26
27 %\renewcommand\citemid{\ } % no comma in optional referenc note
28 \bibpunct{}{}{,}{s}{}{;}
29 \bibliographystyle{aip}
30
31
32 % \documentclass[journal = jpccck, manuscript = article]{achemso}
33 % \setkeys{acs}{usetitle = true}
34 % \usepackage{achemso}
35 % \usepackage{natbib}
36 % \usepackage{multirow}
37 % \usepackage{wrapfig}
38 % \usepackage{fixltx2e}
39
40 \usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions
41 \usepackage{url}
42
43
44 \begin{document}
45
46 \title{Simulations of solid-liquid friction at ice-I$_\mathrm{h}$ /
47 water interfaces}
48
49 \author{Patrick B. Louden and J. Daniel
50 Gezelter\footnote{Corresponding author. \ Electronic mail:
51 gezelter@nd.edu} \\
52 Department of Chemistry and Biochemistry,\\
53 University of Notre Dame\\
54 Notre Dame, Indiana 46556}
55
56 \date{\today}
57 \maketitle
58 \begin{doublespace}
59
60
61 \begin{abstract}
62 We have investigated the structural and dynamic properties of the
63 basal and prismatic facets of the ice I$_\mathrm{h}$ / water
64 interface when the solid phase is drawn through the liquid
65 (i.e. sheared relative to the fluid phase). To impose the shear, we
66 utilized a velocity-shearing and scaling (VSS) approach to reverse
67 non-equilibrium molecular dynamics (RNEMD). This method can create
68 simultaneous temperature and velocity gradients and allow the
69 measurement of transport properties at interfaces. The interfacial
70 width was found to be independent of the relative velocity of the
71 ice and liquid layers over a wide range of shear rates. Decays of
72 molecular orientational time correlation functions gave similar
73 estimates for the width of the interfaces, although the short- and
74 longer-time decay components behave differently closer to the
75 interface. Although both facets of ice are in ``stick'' boundary
76 conditions in liquid water, the solid-liquid friction coefficients
77 were found to be significantly different for the basal and prismatic
78 facets of ice.
79 \end{abstract}
80
81 \newpage
82
83 \section{Introduction}
84 %-----Outline of Intro---------------
85 % in general, ice/water interface is important b/c ....
86 % here are some people who have worked on ice/water, trying to understand the processes above ....
87 % with the recent development of VSS-RNEMD, we can now look at the shearing problem
88 % talk about what we will present in this paper
89 % -------End Intro------------------
90
91 %Gay02: cites many other ice/water papers, make sure to cite them.
92
93 Understanding the ice/water interface is essential for explaining
94 complex processes such as nucleation and crystal
95 growth,\cite{Han92,Granasy95,Vanfleet95} crystal
96 melting,\cite{Weber83,Han92,Sakai96,Sakai96B} and some fascinating
97 biological processes, such as the behavior of the antifreeze proteins
98 found in winter flounder,\cite{Wierzbicki07, Chapsky97} and certain
99 terrestrial arthropods.\cite{Duman:2001qy,Meister29012013} There has
100 been significant progress on understanding the structure and dynamics
101 of quiescent ice/water interfaces utilizing both theory and
102 experiment. Haymet \emph{et al.} have done extensive work on ice I$_\mathrm{h}$,
103 including characterizing and determining the width of the ice/water
104 interface for the SPC,\cite{Karim90} SPC/E,\cite{Gay02,Bryk02} CF1,\cite{Hayward01,Hayward02} and TIP4P~\cite{Karim88} models for
105 water.
106 More recently, Haymet \emph{et al.} have investigated the effects
107 cations and anions have on crystal
108 nucleation.\cite{Bryk04,Smith05,Wilson08,Wilson10} Nada \emph{et al.}
109 have also studied ice/water
110 interfaces,\cite{Nada95,Nada00,Nada03,Nada12} and have found that the
111 differential growth rates of the facets of ice I$_\mathrm{h}$ are due to the the
112 reordering of the hydrogen bonding network.\cite{Nada05}
113
114 The movement of liquid water over the facets of ice has been less
115 thoroughly studied than the quiescent surfaces. This process is
116 potentially important in understanding transport of large blocks of
117 ice in water (which has important implications in the earth sciences),
118 as well as the relative motion of crystal-crystal interfaces that have
119 been separated by nanometer-scale fluid domains. In addition to
120 understanding both the structure and thickness of the interfacial
121 regions, it is important to understand the molecular origin of
122 friction, drag, and other changes in dynamical properties of the
123 liquid in the regions close to the surface that are altered by the
124 presence of a shearing of the bulk fluid relative to the solid phase.
125
126 In this work, we apply a recently-developed velocity shearing and
127 scaling approach to reverse non-equilibrium molecular dynamics
128 (VSS-RNEMD). This method makes it possible to calculate transport
129 properties like the interfacial thermal conductance across
130 heterogeneous interfaces,\cite{Kuang12} and can create simultaneous
131 temperature and velocity gradients and allow the measurement of
132 friction and thermal transport properties at interfaces. This has
133 allowed us to investigate the width of the ice/water interface as the
134 ice is sheared through the liquid, while simultaneously imposing a
135 weak thermal gradient to prevent frictional heating of the interface.
136 In the sections that follow, we discuss the methodology for creating
137 and simulating ice/water interfaces under shear and provide results
138 from both structural and dynamical correlation functions. We also
139 show that the solid-liquid interfacial friction coefficient depends
140 sensitively on the details of the surface morphology.
141
142 \section{Methodology}
143
144 \subsection{Stable ice I$_\mathrm{h}$ / water interfaces under shear}
145
146 The structure of ice I$_\mathrm{h}$ is well understood; it
147 crystallizes in a hexagonal space group P$6_3/mmc$, and the hexagonal
148 crystals of ice have two faces that are commonly exposed, the basal
149 face $\{0~0~0~1\}$, which forms the top and bottom of each hexagonal
150 plate, and the prismatic $\{1~0~\bar{1}~0\}$ face which forms the
151 sides of the plate. Other less-common, but still important, faces of
152 ice I$_\mathrm{h}$ are the secondary prism, $\{1~1~\bar{2}~0\}$, and
153 pyramidal, $\{2~0~\bar{2}~1\}$, faces. Ice I$_\mathrm{h}$ is normally
154 proton disordered in bulk crystals, although the surfaces probably
155 have a preference for proton ordering along strips of dangling H-atoms
156 and Oxygen lone pairs.\cite{Buch:2008fk}
157
158 \begin{table}[h]
159 \centering
160 \caption{Mapping between the Miller indices of four facets of ice in
161 the $P6_3/mmc$ crystal system to the orthorhombic $P2_12_12_1$
162 system in reference \bibpunct{}{}{,}{n}{}{,} \protect\citep{Hirsch04}.}
163 \label{tab:equiv}
164 \begin{tabular}{|ccc|} \hline
165 & hexagonal & orthorhombic \\
166 & ($P6_3/mmc$) & ($P2_12_12_1$) \\
167 crystal face & Miller indices & equivalent \\ \hline
168 basal & $\{0~0~0~1\}$ & $\{0~0~1\}$ \\
169 prism & $\{1~0~\bar{1}~0\}$ & $\{1~0~0\}$ \\
170 secondary prism & $\{1~1~\bar{2}~0\}$ & $\{1~3~0\}$ \\
171 pyramidal & $\{2~0~\bar{2}~1\}$ & $\{2~0~1\}$ \\ \hline
172 \end{tabular}
173 \end{table}
174 For small simulated ice interfaces, it is useful to work with
175 proton-ordered, but zero-dipole crystal that exposes these strips of
176 dangling H-atoms and lone pairs. When placing another material in
177 contact with one of the ice crystalline planes, it is also quite
178 useful to have an orthorhombic (rectangular) box. Recent work by
179 Hirsch and Ojam\"{a}e describes a number of alternative crystal
180 systems for proton-ordered bulk ice I$_\mathrm{h}$ using orthorhombic
181 cells.\cite{Hirsch04}
182
183 In this work, we are using Hirsch and Ojam\"{a}e's structure 6 which
184 is an orthorhombic cell ($P2_12_12_1$) that produces a proton-ordered
185 version of ice I$_\mathrm{h}$. Table \ref{tab:equiv} contains a
186 mapping between the Miller indices of common ice facets in the
187 P$6_3/mmc$ crystal system and those in the Hirsch and Ojam\"{a}e
188 $P2_12_12_1$ system.
189
190 Structure 6 from the Hirsch and Ojam\"{a}e paper has lattice
191 parameters $a = 4.49225$ \AA\ , $b = 7.78080$ \AA\ , $c = 7.33581$ \AA\
192 and two water molecules whose atoms reside at fractional coordinates
193 given in table \ref{tab:p212121}. To construct the basal and prismatic
194 interfaces, these crystallographic coordinates were used to construct
195 an orthorhombic unit cell which was then replicated in all three
196 dimensions yielding a proton-ordered block of ice I$_{h}$. To expose
197 the desired face, the orthorhombic representation was then cut along
198 the ($001$) or ($100$) planes for the basal and prismatic faces
199 respectively. The resulting block was rotated so that the exposed
200 faces were aligned with the $z$-axis normal to the exposed face. The
201 block was then cut along two perpendicular directions in a way that
202 allowed for perfect periodic replication in the $x$ and $y$ axes,
203 creating a slab with either the basal or prismatic faces exposed along
204 the $z$ axis. The slab was then replicated in the $x$ and $y$
205 dimensions until a desired sample size was obtained.
206
207 \begin{table}[h]
208 \centering
209 \caption{Fractional coordinates for water in the orthorhombic
210 $P2_12_12_1$ system for ice I$_\mathrm{h}$ in reference \bibpunct{}{}{,}{n}{}{,} \protect\citep{Hirsch04}.}
211 \label{tab:p212121}
212 \begin{tabular}{|cccc|} \hline
213 atom type & x & y & z \\ \hline
214 O & 0.7500 & 0.1667 & 0.4375 \\
215 H & 0.5735 & 0.2202 & 0.4836 \\
216 H & 0.7420 & 0.0517 & 0.4836 \\
217 O & 0.2500 & 0.6667 & 0.4375 \\
218 H & 0.2580 & 0.6693 & 0.3071 \\
219 H & 0.4265 & 0.7255 & 0.4756 \\ \hline
220 \end{tabular}
221 \end{table}
222
223 Our ice / water interfaces were created using a box of liquid water
224 that had the same dimensions (in $x$ and $y$) as the ice block.
225 Although the experimental solid/liquid coexistence temperature under
226 atmospheric pressure is close to 273~K, Haymet \emph{et al.} have done
227 extensive work on characterizing the ice/water interface, and find
228 that the coexistence temperature for simulated water is often quite a
229 bit different.\cite{Karim88,Karim90,Hayward01,Bryk02,Hayward02} They
230 have found that for the SPC/E water model,\cite{Berendsen87} which is
231 also used in this study, the ice/water interface is most stable at
232 225~$\pm$5~K.\cite{Bryk02} This liquid box was therefore equilibrated at
233 225~K and 1~atm of pressure in the NPAT ensemble (with the $z$ axis
234 allowed to fluctuate to equilibrate to the correct pressure). The
235 liquid and solid systems were combined by carving out any water
236 molecule from the liquid simulation cell that was within 3~\AA\ of any
237 atom in the ice slab. The resulting basal system was $23.87 \times 35.83
238 \times 98.64$ \AA\ with 900 SPC/E molecules in the ice slab, and 1846 in
239 the liquid phase. Similarly, the prismatic system was $36.12 \times 36.43
240 \times 86.10$ \AA\ with 1000 SPC/E molecules in the ice slab and 2684 in
241 the liquid.
242
243 Molecular translation and orientational restraints were applied in the
244 early stages of equilibration to prevent melting of the ice slab.
245 These restraints were removed during NVT equilibration, well before
246 data collection was carried out.
247
248 \subsection{Shearing ice / water interfaces without bulk melting}
249
250 As a solid is dragged through a liquid, there is frictional heating
251 that will act to melt the interface. To study the behavior of the
252 interface under a shear stress without causing the interface to melt,
253 it is necessary to apply a weak thermal gradient in combination with
254 the momentum gradient. This can be accomplished using the velocity
255 shearing and scaling (VSS) variant of reverse non-equilibrium
256 molecular dynamics (RNEMD), which utilizes a series of simultaneous
257 velocity exchanges between two regions within the simulation
258 cell.\cite{Kuang12} One of these regions is centered within the ice
259 slab, while the other is centrally located in the liquid
260 region. VSS-RNEMD provides a set of conservation constraints for
261 creating either a momentum flux or a thermal flux (or both
262 simultaneously) between the two slabs. Satisfying the constraint
263 equations ensures that the new configurations are sampled from the
264 same NVE ensemble as before the VSS move.
265
266 The VSS moves are applied periodically to scale and shift the particle
267 velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in two slabs ($H$ and
268 $C$) which are separated by half of the simulation box,
269 \begin{displaymath}
270 \begin{array}{rclcl}
271
272 & \underline{\mathrm{shearing}} & &
273 \underline{~~~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~~~} \\
274 \mathbf{v}_i \leftarrow &
275 \mathbf{a}_c\ & + & c\cdot\left(\mathbf{v}_i - \langle\mathbf{v}_c
276 \rangle\right) + \langle\mathbf{v}_c\rangle \\
277 \mathbf{v}_j \leftarrow &
278 \mathbf{a}_h & + & h\cdot\left(\mathbf{v}_j - \langle\mathbf{v}_h
279 \rangle\right) + \langle\mathbf{v}_h\rangle .
280
281 \end{array}
282 \end{displaymath}
283 Here $\langle\mathbf{v}_c\rangle$ and $\langle\mathbf{v}_h\rangle$ are
284 the center of mass velocities in the $C$ and $H$ slabs, respectively.
285 Within the two slabs, particles receive incremental changes or a
286 ``shear'' to their velocities. The amount of shear is governed by the
287 imposed momentum flux, $\mathbf{j}_z(\mathbf{p})$
288 \begin{eqnarray}
289 \mathbf{a}_c & = & - \mathbf{j}_z(\mathbf{p}) \Delta t / M_c \label{vss1}\\
290 \mathbf{a}_h & = & + \mathbf{j}_z(\mathbf{p}) \Delta t / M_h \label{vss2}
291 \end{eqnarray}
292 where $M_{\{c,h\}}$ is the total mass of particles within each of the
293 slabs and $\Delta t$ is the interval between two separate operations.
294
295 To simultaneously impose a thermal flux ($J_z$) between the slabs we
296 use energy conservation constraints,
297 \begin{eqnarray}
298 K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\mathbf{v}_c
299 \rangle^2) + \frac{1}{2}M_c (\langle \mathbf{v}_c \rangle + \mathbf{a}_c)^2 \label{vss3}\\
300 K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\mathbf{v}_h
301 \rangle^2) + \frac{1}{2}M_h (\langle \mathbf{v}_h \rangle +
302 \mathbf{a}_h)^2 \label{vss4}.
303 \label{constraint}
304 \end{eqnarray}
305 Simultaneous solution of these quadratic formulae for the scaling
306 coefficients, $c$ and $h$, will ensure that the simulation samples from
307 the original microcanonical (NVE) ensemble. Here $K_{\{c,h\}}$ is the
308 instantaneous translational kinetic energy of each slab. At each time
309 interval, it is a simple matter to solve for $c$, $h$, $\mathbf{a}_c$,
310 and $\mathbf{a}_h$, subject to the imposed momentum flux,
311 $j_z(\mathbf{p})$, and thermal flux, $J_z$, values. Since the VSS
312 operations do not change the kinetic energy due to orientational
313 degrees of freedom or the potential energy of a system, configurations
314 after the VSS move have exactly the same energy (and linear
315 momentum) as before the move.
316
317 As the simulation progresses, the VSS moves are performed on a regular
318 basis, and the system develops a thermal and/or velocity gradient in
319 response to the applied flux. In a bulk material, it is quite simple
320 to use the slope of the temperature or velocity gradients to obtain
321 either the thermal conductivity or shear viscosity.
322
323 The VSS-RNEMD approach is versatile in that it may be used to
324 implement thermal and shear transport simultaneously. Perturbations
325 of velocities away from the ideal Maxwell-Boltzmann distributions are
326 minimal, as is thermal anisotropy. This ability to generate
327 simultaneous thermal and shear fluxes has been previously utilized to
328 map out the shear viscosity of SPC/E water over a wide range of
329 temperatures (90~K) with a single 1~ns simulation.\cite{Kuang12}
330
331 For this work, we are using the VSS-RNEMD method primarily to generate
332 a shear between the ice slab and the liquid phase, while using a weak
333 thermal gradient to maintain the interface at the 225~K target
334 value. This ensures minimal melting of the bulk ice phase and allows
335 us to control the exact temperature of the interface.
336
337 \subsection{Computational Details}
338 All simulations were performed using OpenMD,\cite{OOPSE,openmd} with a
339 time step of 2 fs and periodic boundary conditions in all three
340 dimensions. Electrostatics were handled using the damped-shifted
341 force real-space electrostatic kernel.\cite{Ewald} The systems were
342 divided into 100 bins along the $z$-axis for the VSS-RNEMD moves,
343 which were attempted every 50~fs.
344
345 The interfaces were equilibrated for a total of 10 ns at equilibrium
346 conditions before being exposed to either a shear or thermal gradient.
347 This consisted of 5 ns under a constant temperature (NVT) integrator
348 set to 225K followed by 5 ns under a microcanonical integrator. Weak
349 thermal gradients were allowed to develop using the VSS-RNEMD (NVE)
350 integrator using a small thermal flux ($-2.0\times 10^{-6}$
351 kcal/mol/\AA$^2$/fs) for a duration of 5 ns to allow the gradient to
352 stabilize. The resulting temperature gradient was $\approx$ 10K over
353 the entire 100 \AA\ box length, which was sufficient to keep the
354 temperature at the interface within $\pm 1$ K of the 225K target.
355
356 Velocity gradients were then imposed using the VSS-RNEMD (NVE)
357 integrator with a range of momentum fluxes. These gradients were
358 allowed to stabilize for 1~ns before data collection began. Once
359 established, four successive 0.5~ns runs were performed for each shear
360 rate. During these simulations, snapshots of the system were taken
361 every 1~ps, and statistics on the structure and dynamics in each bin
362 were accumulated throughout the simulations. Although there was some
363 small variation in the measured interfacial width between succcessive
364 runs, no indication of bulk melting or crystallization was observed.
365 That is, no large scale changes in the positions of the left and right
366 interfaces occurred during the simulations.
367
368 \section{Results and discussion}
369
370 \subsection{Interfacial width}
371 Any order parameter or time correlation function that changes as one
372 crosses an interface from a bulk liquid to a solid can be used to
373 measure the width of the interface. In previous work on the ice/water
374 interface, Haymet {\it et al.}\cite{Bryk02} have utilized structural
375 features (including the density) as well as dynamic properties
376 (including the diffusion constant) to estimate the width of the
377 interfaces for a number of facets of the ice crystals. Because
378 VSS-RNEMD imposes a lateral flow, parameters that depend on
379 translational motion of the molecules (e.g. diffusion) may be
380 artificially skewed by the RNEMD moves. A structural parameter is not
381 influenced by the RNEMD perturbations to the same degree. Here, we
382 have used the local tetrahedral order parameter as described by
383 Kumar\cite{Kumar09} and Errington\cite{Errington01} as our principal
384 measure of the interfacial width. A previous study by Bryk and Haymet
385 also used local tetrahedrality as an order parameter for ice/water
386 interfaces.\cite{Bryk2004b}
387
388 The local tetrahedral order parameter, $q(z)$, is given by
389 \begin{equation}
390 q(z) = \int_0^L \sum_{k=1}^{N} \Bigg(1 -\frac{3}{8}\sum_{i=1}^{3}
391 \sum_{j=i+1}^{4} \bigg(\cos\psi_{ikj}+\frac{1}{3}\bigg)^2\Bigg)
392 \delta(z_{k}-z)\mathrm{d}z \Bigg/ N_z
393 \label{eq:qz}
394 \end{equation}
395 where $\psi_{ikj}$ is the angle formed between the oxygen site on
396 central molecule $k$, and the oxygen sites on two of the four closest
397 molecules, $i$ and $j$. Molecules $i$ and $j$ are further restricted
398 to lie withing the first peak in the pair distribution function for
399 molecule $k$ (typically $<$ 3.41 \AA\ for water). $N_z = \int
400 \delta(z_k - z) \mathrm{d}z$ is a normalization factor to account for
401 the varying population of molecules within each finite-width bin. The
402 local tetrahedral order parameter has a range of $(0,1)$, where the
403 larger values of $q$ indicate a larger degree of tetrahedral ordering
404 of the local environment. In perfect ice I$_\mathrm{h}$ structures,
405 the parameter can approach 1 at low temperatures, while in liquid
406 water, the ordering is significantly less tetrahedral, and values of
407 $q(z) \approx 0.75$ are more common.
408
409 To estimate the interfacial width, the system was divided into 100
410 bins along the $z$-dimension. The $q_{z}$ function was time-averaged
411 to give yield a tetrahedrality profile of the system. The profile was
412 then fit to a hyperbolic tangent that smoothly links the liquid and
413 solid states,
414 \begin{equation}\label{tet_fit}
415 q(z) \approx
416 q_{liq}+\frac{q_{ice}-q_{liq}}{2}\left[\tanh\left(\frac{z-l}{w}\right)-\tanh\left(\frac{z-r}{w}\right)\right]+\beta\left|z-
417 \frac{r+l}{2}\right|.
418 \end{equation}
419 Here $q_{liq}$ and $q_{ice}$ are the local tetrahedral order parameter
420 for the bulk liquid and ice domains, respectively, $w$ is the width of
421 the interface. $l$ and $r$ are the midpoints of the left and right
422 interfaces, respectively. The last term in eq. \eqref{tet_fit}
423 accounts for the influence that the weak thermal gradient has on the
424 tetrahedrality profile in the liquid region.
425
426 In Figures \ref{fig:bComic} and \ref{fig:pComic} we see the
427 $z$-coordinate profiles for tetrahedrality, temperature, and the
428 $x$-component of the velocity for the basal and prismatic interfaces.
429 The lower panels show the $q(z)$ (black circles) along with the
430 hyperbolic tangent fits (red lines). In the liquid region, the local
431 tetrahedral order parameter, $q(z) \approx 0.75$ while in the
432 crystalline region, $q(z) \approx 0.94$, indicating a more tetrahedral
433 environment. The vertical dotted lines denote the midpoint of the
434 interfaces ($r$ and $l$ in eq. \eqref{tet_fit}). The weak thermal
435 gradient applied to the systems in order to keep the interface at
436 225~$\pm$~5K, can be seen in middle panels. The transverse velocity
437 profile is shown in the upper panels. It is clear from the upper
438 panels that water molecules in close proximity to the surface (i.e.
439 within 10~\AA\ to 15~\AA~of the interfaces) have transverse velocities
440 quite close to the velocities within the ice block. There is no
441 velocity discontinuity at the interface, which indicates that the
442 shearing of ice/water interfaces occurs in the ``stick'' or no-slip
443 boundary conditions.
444
445 \begin{figure}
446 \includegraphics[width=\linewidth]{bComicStrip}
447 \caption{\label{fig:bComic} The basal interface with a shear rate of
448 1.3 ms\textsuperscript{-1}. Lower panel: the local tetrahedral
449 order parameter, $q(z)$, (black circles) and the hyperbolic tangent
450 fit (red line). Middle panel: the imposed thermal gradient required
451 to maintain a fixed interfacial temperature. Upper panel: the
452 transverse velocity gradient that develops in response to an imposed
453 momentum flux. The vertical dotted lines indicate the locations of
454 the midpoints of the two interfaces.}
455 \end{figure}
456
457 \begin{figure}
458 \includegraphics[width=\linewidth]{pComicStrip}
459 \caption{\label{fig:pComic} The prismatic interface with a shear rate
460 of 2.0 ms\textsuperscript{-1}. Panel
461 descriptions match those in figure \ref{fig:bComic}}
462 \end{figure}
463
464 From the fits using eq. \eqref{tet_fit}, we find the interfacial width
465 for the basal and prismatic systems to be 3.2~$\pm$~0.4~\AA\ and
466 3.6~$\pm$~0.2~\AA\ , respectively, with no applied momentum flux. Over
467 the range of shear rates investigated, $0.6 \pm 0.3 \mathrm{ms}^{-1}
468 \rightarrow 5.3 \pm 0.5 \mathrm{ms}^{-1}$ for the basal system and
469 $0.9 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 4.5 \pm 0.1
470 \mathrm{ms}^{-1}$ for the prismatic, we found no appreciable change in
471 the interface width. The fit values for the interfacial width ($w$)
472 over all shear rates contained the values reported above within their
473 error bars. Note that the interfacial widths reported here are based
474 on the hyperbolic tangent parameter $w$ in Eq. \ref{tet_fit}. This is
475 related to, but not identical with, the 10\%-90\% intefacial widths
476 commonly used in previous studies.\cite{Bryk02,Bryk2004b} To estimate
477 the 10\%-90\% widths, it is a simple matter to scale the widths
478 obtained from the hyperbolic tangent fits to obtain $w_{10-90} =
479 2.197 \times w$.\cite{Bryk02,Bryk2004b} This results in $w_{10-90}$
480 values of 7.0~$\pm$~0.9~\AA\ for the basal face, and 7.9~$\pm$~0.4
481 \AA\ for the prismatic face. These are somewhat smaller than
482 previously reported values.
483
484 \begin{figure}
485 \includegraphics[width=\linewidth]{interface_width_by_shear_rate}
486 \caption{\label{fig:widthByShear} The width of the ice water
487 interfaces (as measured by Eq. \ref{tet_fit}) exhibits no dependence
488 on the applied shear rate between the ice and water regions.}
489 \end{figure}
490
491
492
493 \subsubsection{Orientational Dynamics}
494 The orientational time correlation function,
495 \begin{equation}\label{C(t)1}
496 C_{2}(t)=\langle P_{2}(\mathbf{u}(0) \cdot \mathbf{u}(t)) \rangle,
497 \end{equation}
498 gives insight into the local dynamic environment around the water
499 molecules. The rate at which the function decays provides information
500 about hindered motions and the timescales for relaxation. In
501 eq. \eqref{C(t)1}, $P_{2}$ is the second-order Legendre polynomial,
502 the vector $\mathbf{u}$ is often taken as HOH bisector, although
503 slightly different behavior can be observed when $\mathbf{u}$ is the
504 vector along one of the OH bonds. The angle brackets denote an
505 ensemble average over all water molecules in a given spatial region.
506
507 To investigate the dynamic behavior of water at the ice interfaces, we
508 have computed $C_{2}(z,t)$ for molecules that are present within a
509 particular slab along the $z$- axis at the initial time. The change
510 in the decay behavior as a function of the $z$ coordinate is another
511 measure of the change of how the local environment changes across the
512 ice/water interface. To compute these correlation functions, each of
513 the 0.5 ns simulations was followed by a shorter 200 ps microcanonical
514 (NVE) simulation in which the positions and orientations of every
515 molecule in the system were recorded every 0.1 ps. The systems were
516 then divided into 30 bins along the $z$-axis and $C_2(t)$ was
517 evaluated for each bin.
518
519 In simulations of water at biological interfaces, Furse {\em et al.}
520 fit $C_2(t)$ functions for water with triexponential
521 functions,\cite{Furse08} where the three components of the decay
522 correspond to a fast ($<$200 fs) reorientational piece driven by the
523 restoring forces of existing hydrogen bonds, a middle (on the order of
524 several ps) piece describing the large angle jumps that occur during
525 the breaking and formation of new hydrogen bonds,and a slow (on the
526 order of tens of ps) contribution describing the translational motion
527 of the molecules. The model for orientational decay presented
528 recently by Laage and Hynes {\em et al.}\cite{Laage08,Laage11} also
529 includes three similar decay constants, although two of the time
530 constants are linked, and the resulting decay curve has two parameters
531 governing the dynamics of decay.
532
533 In our ice/water interfaces, we are at substantially lower
534 temperatures, and the water molecules are further perturbed by the
535 presence of the ice phase nearby. We have obtained the most
536 reasonable fits using triexponential functions with three distinct
537 time domains, as well as a constant piece to account for the water
538 stuck in the ice phase that does not experience any long-time
539 orientational decay,
540 \begin{equation}
541 C_{2}(t) \approx a e^{-t/\tau_\mathrm{short}} + b e^{-t/\tau_\mathrm{middle}} + c
542 e^{-t/\tau_\mathrm{long}} + (1-a-b-c)
543 \end{equation}
544 Average values for the three decay constants (and error estimates)
545 were obtained for each bin. In figures \ref{fig:basal_Tau_comic_strip}
546 and \ref{fig:prismatic_Tau_comic_strip}, the three orientational decay
547 times are shown as a function of distance from the center of the ice
548 slab.
549
550 \begin{figure}
551 \includegraphics[width=\linewidth]{basal_Tau_comic_strip}
552 \caption{\label{fig:basal_Tau_comic_strip} The three decay constants
553 of the orientational time correlation function, $C_2(t)$, for water
554 as a function of distance from the center of the ice slab. The
555 dashed line indicates the location of the basal face (as determined
556 from the tetrahedrality order parameter) and the black and red lines
557 are fits using Eq. \ref{tauFit}. The moderate and long time
558 contributions slow down close to the interface which would be
559 expected under reorganizations that involve large motions of the
560 molecules (e.g. frame-reorientations and jumps). The observed
561 speed-up in the short time contribution is surprising, but appears
562 to reflect the restricted motion of librations closer to the
563 interface.}
564 \end{figure}
565
566 \begin{figure}
567 \includegraphics[width=\linewidth]{prismatic_Tau_comic_strip}
568 \caption{\label{fig:prismatic_Tau_comic_strip}
569 Decay constants for $C_2(t)$ at the prismatic interface. Panel
570 descriptions match those in figure \ref{fig:basal_Tau_comic_strip}.}
571 \end{figure}
572
573 Figures \ref{fig:basal_Tau_comic_strip} and
574 \ref{fig:prismatic_Tau_comic_strip} show the three decay constants for
575 the orientational time correlation function for water at varying
576 displacements from the center of the ice slab for both the basal and
577 prismatic interfaces. The vertical dotted lines indicate the
578 locations of the midpoints of the interfaces as determined by the
579 tetrahedrality fits. In the liquid regions, $\tau_\mathrm{middle}$ and
580 $\tau_\mathrm{long}$ have consistent values around 3-4 ps and 20-40 ps,
581 respectively, and increase in value approaching the interface.
582 According to the jump model of Laage and Hynes {\em et
583 al.},\cite{Laage08,Laage11} $\tau_\mathrm{middle}$ corresponds to the
584 breaking and making of hydrogen bonds and $\tau_\mathrm{long}$ is explained
585 with translational motion of the molecules (i.e. frame reorientation).
586 The shortest of the three decay constants, the librational time
587 $\tau_\mathrm{short}$ has a value of about 70 fs in the liquid region,
588 and decreases in value approaching the interface. The observed
589 speed-up in the short time contribution is surprising, but appears to
590 reflect the restricted motion of librations closer to the interface.
591
592 The control systems (with no applied momentum flux) are shown with
593 black symbols in figs. \ref{fig:basal_Tau_comic_strip} and
594 \ref{fig:prismatic_Tau_comic_strip}, while those obtained while a
595 shear was active are shown in red. There are two notable features
596 from the spatial dependence of the orientational correlation times
597 seen in figures \ref{fig:basal_Tau_comic_strip} and
598 \ref{fig:prismatic_Tau_comic_strip}. First, there are nearly-constant
599 liquid-state values for $\tau_\mathrm{short}$, $\tau_\mathrm{middle}$,
600 and $\tau_\mathrm{long}$ at large displacements from the
601 interface. Second, there appears to be a single distance, $d$ from the
602 interface at which all three decay times begin to deviate from their
603 bulk liquid values. To quantify this ``dynamic interface width,''
604 each of the decay constant profiles were fit to
605 \begin{equation}\label{tauFit}
606 \tau_{a}(z) = \tau^{\mathrm{liquid}}_{a}+(\tau^{\mathrm{wall}}_a-\tau^{\mathrm{liquid}}_a)e^{-(z-z^{\mathrm{wall}})/d}
607 \end{equation}
608 where $a = \{ \mathrm{short, middle,} ~or~ \mathrm{long} \}$. In this
609 expression, $\tau^{\mathrm{liquid}}_a$ and $\tau^{\mathrm{wall}}_a$
610 are the bulk liquid and projected values of the decay constants at the
611 wall, $z^{\mathrm{wall}}$ is the location of the interface, and $d$ is
612 the dynamic interface width. The widths for the two facets of ice,
613 $d_{basal}$ and $d_{prismatic}$, were determined by averaging $d$
614 values for each of the three decay times. For the basal system, we
615 find $d_{basal}$ to be 2.9 \AA\ for the quiescent control set, and 2.8
616 \AA\ for a simulation with a shear rate of 1.3
617 ms\textsuperscript{-1}. We also find $d_{prismatic}$ to be slightly
618 larger than $d_{basal}$ for both the control and an applied shear,
619 with widths of 3.6 \AA\ for the quiescent system and 3.5 \AA\ for a
620 simulation with a 2 ms\textsuperscript{-1} shear rate. There appears
621 to be little or no dependence of the dynamic interface width on the
622 applied shear rate.
623
624 The spatial dependence of the decay times in $C_2(t)$ determines a
625 dynamic interfacial width that is in good agreement with the
626 structural interfacial width determined by the local tetrahedrality.
627 For exponential models like the one in Eq. \ref{tauFit}, it is again a
628 simple matter to estimate the 10\%-90\% interfacial width by
629 multiplying the $d$ constant by $2.19723$. For the interfaces
630 studied, we obtain 10\%-90\% dynamic interface width $(d_{10-90})$
631 values of $\sim$~6.3~\AA\ for the basal face and 7.8~\AA\ for the
632 prismatic face.
633
634 \subsection{Coefficient of Friction of the Interface}
635 As liquid water flows over an ice interface, there is a distance from
636 the structural interface where bulk-like hydrodynamics are recovered.
637 Bocquet and Barrat constructed a theory for the hydrodynamic boundary
638 parameters, which include the slipping length
639 $\left(\delta_\mathrm{wall}\right)$ of this boundary layer and the
640 ``hydrodynamic position'' of the boundary
641 $\left(z_\mathrm{wall}\right)$.\cite{PhysRevLett.70.2726,PhysRevE.49.3079}
642 This last parameter is the location (relative to a solid surface)
643 where the bulk-like behavior is recovered. Work by Mundy {\it et al.}
644 has helped to combine these parameters into a liquid-solid friction
645 coefficient, which quantifies the resistance to pulling the solid
646 interface through a liquid,\cite{Mundy1997305}
647 \begin{equation}
648 \lambda_\mathrm{wall} = \frac{\eta}{\delta_\mathrm{wall}}.
649 \end{equation}
650 This expression is nearly identical to one provided by Pit {\it et
651 al.} for the solid-liquid friction of an interface,\cite{Pit99}
652 \begin{equation}\label{Pit}
653 \lambda=\frac{\eta}{\delta}
654 \end{equation}
655 where $\delta$ is the slip length for the liquid measured at the
656 location of the interface itself. In our simulations, the shoulder on
657 the velocity profile indicating the location of the hydrodynamic
658 boundary in the liquid is not always apparent. In some cases, the
659 linear behavior persists nearly up to the interfacial region. For
660 this reason, the hydrodynamic position of the boundary is not always
661 computable, while the Pit approach (Eq. \ref{Pit}) can be used to find
662 the solid-liquid friction coefficient more reliably.
663
664 In both the Pit and hydrodynamic boundary expressions, $\eta$ is the
665 shear viscosity of the bulk-like region of the liquid, a quantity
666 which is easily obtained in VSS-RNEMD simulations by fitting the
667 velocity profile in the region far from the surface.\cite{Kuang12}
668 Assuming linear response in the bulk-like region,
669 \begin{equation}\label{Kuang}
670 j_{z}(p_{x})=-\eta \left(\frac{\partial v_{x}}{\partial z}\right)
671 \end{equation}
672 Substituting this result into eq. \eqref{Pit}, we can estimate the
673 solid-liquid coefficient using the slip length,
674 \begin{equation}
675 \lambda=-\frac{j_{z}(p_{x})} {\left(\frac{\partial v_{x}}{\partial
676 z}\right) \delta}
677 \end{equation}
678
679 For ice / water interfaces, the boundary conditions are no-slip, so
680 projecting the bulk liquid state velocity profile yields a negative
681 slip length. This length is the difference between the structural edge
682 of the ice (determined by the tetrahedrality profile) and the location
683 where the projected velocity of the bulk liquid intersects the solid
684 phase velocity (see Figure \ref{fig:delta_example}). The coefficients
685 of friction for the basal and the prismatic facets were determined for
686 shearing along both the $x$ and $y$ axes. The values are given in
687 table \ref{tab:lambda}.
688
689 Note that the measured friction coefficient for the basal face is
690 twice that of the prismatic face (regardless of drag direction).
691 These results may seem surprising as the basalface appears smoother
692 than the prismatic with only small undulations of the oxygen
693 positions, while the prismatic surface has deep corrugated channels
694 along the $x$ direction in the crystal system used in this work.
695 However, the corrugations are relatively thin, and the liquid phase
696 water does not appear to populate the channels. The prismatic face
697 therefore effectively presents stripes of solid-phase molecules
698 (making up approximately half of the exposed surface area) with nearly
699 empty space between them. The interfacial friction appears to be
700 independent of the drag direction, so flow parallel to these channels
701 does not explain the lower friction of the prismatic face. A more
702 likely explanation is that the effective contact between the liquid
703 phase and the prismatic face is reduced by the empty corrugations.
704
705 \begin{figure}
706 \includegraphics[width=\linewidth]{delta_example}
707 \caption{\label{fig:delta_example} Determining the (negative) slip
708 length ($\delta$) for the ice-water interfaces (which have decidedly
709 non-slip behavior). This length is the difference between the
710 structural edge of the ice (determined by the tetrahedrality
711 profile) and the location where the projected velocity of the bulk
712 liquid (dashed red line) intersects the solid phase velocity (solid
713 black line). The dotted line indicates the location of the ice as
714 determined by the tetrahedrality profile. This example is taken
715 from the basal-face simulation with an applied shear rate of 3.0
716 ms\textsuperscript{-1}.}
717 \end{figure}
718
719
720 \begin{table}[h]
721 \centering
722 \caption{Solid-liquid friction coefficients (measured in amu~fs\textsuperscript{-1}) }
723 \label{tab:lambda}
724 \begin{tabular}{|ccc|} \hline
725 & \multicolumn{2}{c|}{Drag direction} \\
726 Interface & $x$ & $y$ \\ \hline
727 basal & $0.08 \pm 0.02$ & $0.09 \pm 0.03$ \\
728 prismatic & $0.037 \pm 0.008$ & $0.04 \pm 0.01$ \\ \hline
729 \end{tabular}
730 \end{table}
731
732
733 \section{Conclusion}
734 We have simulated the basal and prismatic facets of an SPC/E model of
735 the ice I$_\mathrm{h}$ / water interface. Using VSS-RNEMD, the ice
736 was sheared relative to the liquid while simultaneously being exposed
737 to a weak thermal gradient which kept the interface at a stable
738 temperature. Calculation of the local tetrahedrality order parameter
739 has shown an apparent independence of the interfacial width on the
740 shear rate. This width was found to be 3.2~$\pm$0.4~\AA\ and
741 3.6~$\pm$0.2~\AA\ for the basal and prismatic systems, respectively.
742 Using the 10\%-90\% width estimates common in previous studies, these
743 interfacial widths are 7.0~$\pm$~0.9~\AA\ for the basal face, and
744 7.9~$\pm$~0.4 \AA\ for the prismatic face.
745
746 Orientational time correlation functions were calculated at varying
747 displacements from the interface, and were found to be similarly
748 independent of the applied momentum flux. The short decay due to the
749 restoring forces of existing hydrogen bonds decreased close to the
750 interface, while the longer-time decay constants increased in close
751 proximity to the interface. There is also an apparent dynamic
752 interface width, $d_{basal}$ and $d_{prismatic}$, at which these
753 deviations from bulk liquid values begin. We found $d_{basal}$ and
754 $d_{prismatic}$ to be approximately 2.8~\AA\ and 3.5~\AA\ . Using the
755 10\%-90\% estimates common in previous simulation studies, the dynamic
756 interface widths are $\sim$~6.3~\AA\ for the basal face and 7.8~\AA\
757 for the prismatic face. These interfacial widths are in good
758 agreement with values determined by the structural analysis of the
759 interface using the local tetrahedrality order parameter.
760
761 The coefficient of liquid-solid friction for each of the facets was
762 also determined and were found to be essentially independent of shear
763 direction relative to the surface orientation. We attribute the
764 factor of two difference between the basal and prismatic friction
765 coefficients to the corrugation of the prismatic face and the
766 resulting reduction of effective surface area for interaction with the
767 liquid.
768
769 \section{Acknowledgements}
770 The authors thank the reviewers for helpful comments and suggestions.
771 Support for this project was provided by the National Science
772 Foundation under grant CHE-0848243. Computational time was provided by
773 the Center for Research Computing (CRC) at the University of Notre
774 Dame.
775
776 \newpage
777 \bibliography{iceWater}
778
779 \end{doublespace}
780
781 % \begin{tocentry}
782 % \begin{wrapfigure}{l}{0.5\textwidth}
783 % \begin{center}
784 % \includegraphics[width=\linewidth]{SystemImage.png}
785 % \end{center}
786 % \end{wrapfigure}
787 % The cell used to simulate liquid-solid shear in ice I$_\mathrm{h}$ /
788 % water interfaces. Velocity gradients were applied using the velocity
789 % shearing and scaling variant of reverse non-equilibrium molecular
790 % dynamics (VSS-RNEMD) with a weak thermal gradient to prevent melting.
791 % The interface width is relatively robust in both structual and dynamic
792 % measures as a function of the applied shear.
793 % \end{tocentry}
794
795 \end{document}
796
797 %**************************************************************
798 %Non-mass weighted slopes (amu*Angstroms^-2 * fs^-1)
799 % basal: slope=0.090677616, error in slope = 0.003691743
800 %prismatic: slope = 0.050101506, error in slope = 0.001348181
801 %Mass weighted slopes (Angstroms^-2 * fs^-1)
802 %basal slope = 4.76598E-06, error in slope = 1.94037E-07
803 %prismatic slope = 3.23131E-06, error in slope = 8.69514E-08
804 %**************************************************************