1 |
% ****** Start of file aipsamp.tex ****** |
2 |
% |
3 |
% This file is part of the AIP files in the AIP distribution for REVTeX 4. |
4 |
% Version 4.1 of REVTeX, October 2009 |
5 |
% |
6 |
% Copyright (c) 2009 American Institute of Physics. |
7 |
% |
8 |
% See the AIP README file for restrictions and more information. |
9 |
% |
10 |
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed |
11 |
% as well as the rest of the prerequisites for REVTeX 4.1 |
12 |
% |
13 |
% It also requires running BibTeX. The commands are as follows: |
14 |
% |
15 |
% 1) latex aipsamp |
16 |
% 2) bibtex aipsamp |
17 |
% 3) latex aipsamp |
18 |
% 4) latex aipsamp |
19 |
% |
20 |
% Use this file as a source of example code for your aip document. |
21 |
% Use the file aiptemplate.tex as a template for your document. |
22 |
\documentclass[% |
23 |
aip,jcp, |
24 |
amsmath,amssymb, |
25 |
preprint,% |
26 |
% reprint,% |
27 |
%author-year,% |
28 |
%author-numerical,% |
29 |
jcp]{revtex4-1} |
30 |
|
31 |
\usepackage{graphicx}% Include figure files |
32 |
\usepackage{dcolumn}% Align table columns on decimal point |
33 |
%\usepackage{bm}% bold math |
34 |
\usepackage{times} |
35 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
36 |
\usepackage{url} |
37 |
|
38 |
\begin{document} |
39 |
|
40 |
\title{Friction at Water / Ice-I$_\mathrm{h}$ interfaces: Do the |
41 |
Facets of Ice Have Different Hydrophilicity?} |
42 |
|
43 |
\author{Patrick B. Louden} |
44 |
|
45 |
\author{J. Daniel Gezelter} |
46 |
\email{gezelter@nd.edu.} |
47 |
\affiliation{Department of Chemistry and Biochemistry, University |
48 |
of Notre Dame, Notre Dame, IN 46556} |
49 |
|
50 |
\date{\today} |
51 |
|
52 |
\begin{abstract} |
53 |
In this follow up paper of the basal and prismatic facets of the |
54 |
Ice-I$_\mathrm{h}$/water interface, we present the |
55 |
pyramidal and secondary prismatic |
56 |
interfaces for both the quiescent and sheared systems. The structural and |
57 |
dynamic interfacial widths for all four crystal facets were found to be in good |
58 |
agreement, and were found to be independent of the shear rate over the shear |
59 |
rates investigated. |
60 |
Decomposition of the molecular orientational time correlation function showed |
61 |
different behavior for the short- and longer-time decay components approaching |
62 |
normal to the interface. Lastly we show through calculation of the interfacial |
63 |
friction coefficient that the basal and pyramidal facets are more |
64 |
hydrophilic than the prismatic and secondary prismatic facets. |
65 |
|
66 |
\end{abstract} |
67 |
|
68 |
\maketitle |
69 |
|
70 |
\section{Introduction} |
71 |
|
72 |
% Why are ice and ice/water interfaces important? |
73 |
There is significant interest in the properties of ice/ice and ice/water |
74 |
interfaces in the geophysics community. Most commonly, the results of shearing |
75 |
two ice blocks past one |
76 |
another\cite{Casassa91, Sukhorukov13, Pritchard12, Lishman13} or the shearing |
77 |
of ice through water\cite{Cuffey99, Bell08}. Using molecular dynamics |
78 |
simulations, Samadashvili has recently shown that when two smooth ice slabs |
79 |
slide past one another, a stable liquid-like layer develops between them.To |
80 |
fundamentally understand these processes an understanding of the ice/water |
81 |
interfaces is needed. |
82 |
|
83 |
Investigation of the ice/water interface is also crucial in understanding |
84 |
the fundamental processes of nucleation, crystal |
85 |
growth,\cite{Han92, Granasy95, Vanfleet95} and crystal |
86 |
melting\cite{Weber83, Han92, Sakai96, Sakai96B}. Insight gained to these |
87 |
properties can also be applied to biological systems of interest, such as |
88 |
the behavior of the antifreeze protein found in winter |
89 |
flounder,\cite{Wierzbicki07,Chapsky97} and certain terrestial |
90 |
arthropods.\cite{Duman:2001qy,Meister29012013} Elucidating the properties which |
91 |
give rise to these processes through experimental techniques can be expensive, |
92 |
complicated, and sometimes infeasible. However, through the use of molecular |
93 |
dynamics simulations much of the problems of investigating these properties |
94 |
are alleviated. |
95 |
|
96 |
Understanding ice/water interfaces inherently begins with the isolated |
97 |
systems. There has been extensive work parameterizing models for liquid water, |
98 |
such as the SPC\cite{Berendsen81}, SPC/E\cite{Berendsen87}, |
99 |
TIP4P\cite{Jorgensen85}, TIP4P/2005\cite{Abascal05}, |
100 |
($\dots$), and more recently to improve these models for simulating |
101 |
the solid phases of water, such as the TIP4P/Ice\cite{Abascal05b} model. The |
102 |
melting point of various crystal structures of ice have been calculated for |
103 |
many of these models |
104 |
(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}), |
105 |
and the partial or complete phase diagram for the model has been determined |
106 |
(SPC/E\cite{Baeaz95,Bryk04b,Sanz04b}, TIP4P\cite{Sanz04,Sanz04b,Koyama04}, TIP5P\cite{Sanz04,Koyama04}). |
107 |
Knowing the melting point for these models has allowed for the investigation |
108 |
of ice/water interfaces. |
109 |
|
110 |
The Ice-I$_\mathrm{h}$/water quiescent interface has been extensively studied |
111 |
over the past 30 years by theory and experiment. Haymet \emph{et al.} have |
112 |
done significant work characterizing and quantifying the width of these |
113 |
interfaces for the SPC,\cite{Karim90} SPC/E,\cite{Gay02,Bryk02}, |
114 |
CF1,\cite{Hayward01,Hayward02} and TIP4P\cite{Karim88} models for water. In |
115 |
recent years, Haymet has focused on investigating the effects cations and |
116 |
anions have on crystal nucleaion and |
117 |
melting.\cite{Bryk04,Smith05,Wilson08,Wilson10} Nada and Furukawa have studied |
118 |
the the basal- and prismatic-water interface width\cite{Nada95}, surface |
119 |
restructuring at temperatures approaching the melting point\cite{Nada00}, |
120 |
and the mechanism of ice growth inhibition by antifreeze |
121 |
proteins\cite{Nada08,Nada11,Nada12}. Nada has developed a six-site water model |
122 |
for ice/water interfaces near the melting point\cite{Nada03}, and studied the |
123 |
dependence of crystal growth shape on applied pressure\cite{Nada11b}. Using |
124 |
this model, Nada and Furukawa have established differential |
125 |
growth rates for the basal, prismatic, and secondary prismatic facets of |
126 |
Ice-I$_\mathrm{h}$ and found they were due to a reordering of the hydrogen |
127 |
bond network in water near the interface\cite{Nada05}. While the work |
128 |
described so far has mainly focused on bulk water on ice, there is significant |
129 |
interest in thin films of water on ice surfaces as well. |
130 |
|
131 |
It is well known that the surface of ice exhibits a premelting layer at |
132 |
temperatures near the melting point, often called a quasi-liquid layer (QLL). |
133 |
Molecular dynamics simulations of the facets of ice-I$_\mathrm{h}$ exposed |
134 |
to vacuum performed by Conde, Vega and Patrykiejew have found QLL widths of |
135 |
approximately 10 \AA\ at 3 K below the melting point\cite{Conde08}. |
136 |
Similarly, Limmer and Chandler have used course grain simulations and |
137 |
statistical field theory to estimated QLL widths at the same temperature to |
138 |
be about 3 nm\cite{Limmer14}. |
139 |
Recently, Sazaki and Furukawa have developed an experimental technique with |
140 |
sufficient spatial and temporal resolution to visulaize and quantitatively |
141 |
analyze QLLs on ice crystals at temperatures near melting\cite{Sazaki10}. They |
142 |
have found the width of the QLLs perpindicular to the surface at -2.2$^{o}$C |
143 |
to be $\mathcal{O}$(\AA). They have also seen the formation of two immiscible |
144 |
QLLs, which displayed different stabilities and dynamics on the crystal |
145 |
surface\cite{Sazaki12}. Knowledge of the hydrophilicities of each |
146 |
of the crystal facets would help further our understanding of the properties |
147 |
and dynamics of the QLLs. |
148 |
|
149 |
Presented here is the follow up to our previous paper\cite{Louden13}, in which |
150 |
the basal and prismatic facets of an ice-I$_\mathrm{h}$/water interface were |
151 |
investigated where the ice was sheared relative to the liquid. By using a |
152 |
recently developed velocity shearing and scaling approach to reverse |
153 |
non-equilibrium molecular dynamics (VSS-RNEMD), simultaneous temperature and |
154 |
velocity gradients can be applied to the system, which allows for measurment |
155 |
of friction and thermal transport properties while maintaining a stable |
156 |
interfacial temperature\cite{Kuang12}. Structural analysis and dynamic |
157 |
correlation functions were used to probe the interfacial response to a shear, |
158 |
and the resulting solid/liquid kinetic friction coefficients were reported. |
159 |
In this paper we present the same analysis for the pyramidal and secondary |
160 |
prismatic facets, and show that the differential interfacial friction |
161 |
coefficients for the four facets of ice-I$_\mathrm{h}$ are determined by their |
162 |
relative hydrophilicity by means of dynamics water contact angle simulations. |
163 |
|
164 |
\section{Methodology} |
165 |
|
166 |
\begin{figure} |
167 |
\includegraphics[width=\linewidth]{Pyr_comic_strip} |
168 |
\caption{\label{fig:pyrComic} The pyramidal interface with a shear |
169 |
rate of 3.8 ms\textsuperscript{-1}. Lower panel: the local tetrahedral order |
170 |
parameter, $q(z)$, (black circles) and the hyperbolic tangent fit (red line). |
171 |
Middle panel: the imposed thermal gradient required to maintain a fixed |
172 |
interfacial temperature. Upper panel: the transverse velocity gradient that |
173 |
develops in response to an imposed momentum flux. The vertical dotted lines |
174 |
indicate the locations of the midpoints of the two interfaces.} |
175 |
\end{figure} |
176 |
|
177 |
\begin{figure} |
178 |
\includegraphics[width=\linewidth]{SP_comic_strip} |
179 |
\caption{\label{fig:spComic} The secondary prismatic interface with a shear |
180 |
rate of 3.5 \ |
181 |
ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:pyrComic}.} |
182 |
\end{figure} |
183 |
|
184 |
\subsection{Pyramidal and secondary prismatic system construction} |
185 |
|
186 |
The construction of the pyramidal and secondary prismatic systems follows that |
187 |
of |
188 |
the basal and prismatic systems presented elsewhere\cite{Louden13}, however |
189 |
the ice crystals and water boxes were equilibrated and combined at 50K |
190 |
instead of 225K. The ice / water systems generated were then equilibrated |
191 |
to 225K. The resulting pyramidal system was |
192 |
$37.47 \times 29.50 \times 93.02$ \AA\ with 1216 |
193 |
SPC/E\cite{Berendsen97} molecules in the ice slab, and 2203 in the liquid |
194 |
phase. The secondary |
195 |
prismatic system generated was $71.87 \times 31.66 \times 161.55$ \AA\ with |
196 |
3840 |
197 |
SPC/E molecules in the ice slab and 8176 molecules in the liquid phase. |
198 |
|
199 |
\subsection{Computational details} |
200 |
% Do we need to justify the sims at 225K? |
201 |
% No crystal growth or shrinkage over 2 successive 1 ns NVT simulations for |
202 |
% either the pyramidal or sec. prismatic ice/water systems. |
203 |
|
204 |
The computational details performed here were equivalent to those reported |
205 |
in our previous publication\cite{Louden13}. The only changes made to the |
206 |
previously reported procedure were the following. VSS-RNEMD moves were |
207 |
attempted every 2 fs instead of every 50 fs. This was done to minimize |
208 |
the magnitude of each individual VSS-RNEMD perturbation to the system. |
209 |
|
210 |
All pyramidal simulations were performed under the canonical (NVT) ensamble |
211 |
except those |
212 |
during which statistics were accumulated for the orientational correlation |
213 |
function, which were performed under the microcanonical (NVE) ensamble. All |
214 |
secondary prismatic |
215 |
simulations were performed under the NVE ensamble. |
216 |
|
217 |
\section{Results and discussion} |
218 |
\subsection{Interfacial width} |
219 |
In the literature there is good agreement that between the solid ice and |
220 |
the bulk water, there exists a region of 'slush-like' water molecules. |
221 |
In this region, the water molecules are structurely distinguishable and |
222 |
behave differently than those of the solid ice or the bulk water. |
223 |
The characteristics of this region have been defined by both structural |
224 |
and dynamic properties; and its width has been measured by the change of these |
225 |
properties from their bulk liquid values to those of the solid ice. |
226 |
Examples of these properties include the density, the diffusion constant, and |
227 |
the translational order profile. \cite{Bryk02,Karim90,Gay02,Hayward01,Hayward02,Karim88} |
228 |
|
229 |
Since the VSS-RNEMD moves used to impose the thermal and velocity gradients |
230 |
perturb the momenta of the water molecules in |
231 |
the systems, parameters that depend on translational motion may give |
232 |
faulty results. A stuructural parameter will be less effected by the |
233 |
VSS-RNEMD perturbations to the system. Due to this, we have used the |
234 |
local tetrahedral order parameter to quantify the width of the interface, |
235 |
which was originally described by Kumar\cite{Kumar09} and |
236 |
Errington\cite{Errington01}, and used by Bryk and Haymet in a previous study |
237 |
of ice/water interfaces.\cite{Bryk2004b} |
238 |
|
239 |
The local tetrahedral order parameter, $q(z)$, is given by |
240 |
\begin{equation} |
241 |
q(z) = \int_0^L \sum_{k=1}^{N} \Bigg(1 -\frac{3}{8}\sum_{i=1}^{3} |
242 |
\sum_{j=i+1}^{4} \bigg(\cos\psi_{ikj}+\frac{1}{3}\bigg)^2\Bigg) |
243 |
\delta(z_{k}-z)\mathrm{d}z \Bigg/ N_z |
244 |
\label{eq:qz} |
245 |
\end{equation} |
246 |
where $\psi_{ikj}$ is the angle formed between the oxygen sites of molecules |
247 |
$i$,$k$, and $j$, where the centeral oxygen is located within molecule $k$ and |
248 |
molecules $i$ and $j$ are two of the closest four water molecules |
249 |
around molecule $k$. All four closest neighbors of molecule $k$ are also |
250 |
required to reside within the first peak of the pair distribution function |
251 |
for molecule $k$ (typically $<$ 3.41 \AA\ for water). |
252 |
$N_z = \int\delta(z_k - z) \mathrm{d}z$ is a normalization factor to account |
253 |
for the varying population of molecules within each finite-width bin. |
254 |
|
255 |
To determine the width of the interfaces, each of the systems were |
256 |
divided into 100 artificial bins along the |
257 |
$z$-dimension, and the local tetrahedral order parameter, $q(z)$, was |
258 |
time-averaged for each of the bins, resulting in a tetrahedrality profile of |
259 |
the system. These profiles are shown across the $z$-dimension of the systems |
260 |
in panel $a$ of Figures \ref{fig:pyrComic} |
261 |
and \ref{fig:spComic} (black circles). The $q(z)$ function has a range of |
262 |
(0,1), where a larger value indicates a more tetrahedral environment. |
263 |
The $q(z)$ for the bulk liquid was found to be $\approx $ 0.77, while values of |
264 |
$\approx $ 0.92 were more common for the ice. The tetrahedrality profiles were |
265 |
fit using a hyperbolic tangent\cite{Louden13} designed to smoothly fit the |
266 |
bulk to ice |
267 |
transition, while accounting for the thermal influence on the profile by the |
268 |
kinetic energy exchanges of the VSS-RNEMD moves. In panels $b$ and $c$, the |
269 |
resulting thermal and velocity gradients from the imposed kinetic energy and |
270 |
momentum fluxes can be seen. The verticle dotted |
271 |
lines traversing all three panels indicate the midpoints of the interface |
272 |
as determined by the hyperbolic tangent fit of the tetrahedrality profiles. |
273 |
|
274 |
From fitting the tetrahedrality profiles for each of the 0.5 nanosecond |
275 |
simulations (panel c of Figures \ref{fig:spComic} and \ref{fig:pyrComic}) |
276 |
by Eq. 6\cite{Louden13},we find the interfacial width to be |
277 |
3.2 $\pm$ 0.2 and 3.2 $\pm$ 0.2 \AA\ for the control system with no applied |
278 |
momentum flux for both the pyramidal and secondary prismatic systems. |
279 |
Over the range of shear rates investigated, |
280 |
0.6 $\pm$ 0.2 $\mathrm{ms}^{-1} \rightarrow$ 5.6 $\pm$ 0.4 $\mathrm{ms}^{-1}$ |
281 |
for the pyramidal system and 0.9 $\pm$ 0.3 $\mathrm{ms}^{-1} \rightarrow$ 5.4 |
282 |
$\pm$ 0.1 $\mathrm{ms}^{-1}$ for the secondary prismatic, we found no |
283 |
significant change in the interfacial width. This follows our previous |
284 |
findings of the basal and |
285 |
prismatic systems, in which the interfacial width was invarient of the |
286 |
shear rate of the ice. The interfacial width of the quiescent basal and |
287 |
prismatic systems was found to be 3.2 $\pm$ 0.4 \AA\ and 3.6 $\pm$ 0.2 \AA\ |
288 |
respectively, over the range of shear rates investigated, 0.6 $\pm$ 0.3 |
289 |
$\mathrm{ms}^{-1} \rightarrow$ 5.3 $\pm$ 0.5 $\mathrm{ms}^{-1}$ for the basal |
290 |
system and 0.9 $\pm$ 0.2 $\mathrm{ms}^{-1} \rightarrow$ 4.5 $\pm$ 0.1 |
291 |
$\mathrm{ms}^{-1}$ for the prismatic. |
292 |
|
293 |
These results indicate that the surface structure of the exposed ice crystal |
294 |
has little to no effect on how far into the bulk the ice-like structural |
295 |
ordering is. Also, it appears that the interface is not structurally effected |
296 |
by shearing the ice through water. |
297 |
|
298 |
|
299 |
\subsection{Orientational dynamics} |
300 |
%Should we include the math here? |
301 |
The orientational time correlation function, |
302 |
\begin{equation}\label{C(t)1} |
303 |
C_{2}(t)=\langle P_{2}(\mathbf{u}(0)\cdot \mathbf{u}(t))\rangle, |
304 |
\end{equation} |
305 |
helps indicate the local environment around the water molecules. The function |
306 |
begins with an initial value of unity, and decays to zero as the water molecule |
307 |
loses memory of its former orientation. Observing the rate at which this decay |
308 |
occurs can provide insight to the mechanism and timescales for the relaxation. |
309 |
In eq. \eqref{C(t)1}, $P_{2}$ is the second-order Legendre polynomial, and |
310 |
$\mathbf{u}$ is the bisecting HOH vector. The angle brackets indicate |
311 |
an ensemble average over all the water molecules in a given spatial region. |
312 |
|
313 |
To investigate the dynamics of the water molecules across the interface, the |
314 |
systems were divided in the $z$-dimension into bins, each $\approx$ 3 \AA\ |
315 |
wide, and \eqref{C(t)1} was computed for each of the bins. A water |
316 |
molecule was allocated to a particular bin if it was initially in the bin |
317 |
at time zero. To compute \eqref{C(t)1}, each 0.5 ns simulation was followed |
318 |
by an additional 200 ps NVE simulation during which the |
319 |
position and orientations of each molecule were recorded every 0.1 ps. |
320 |
|
321 |
The data obtained for each bin was then fit to a triexponential decay given by |
322 |
\begin{equation}\label{C(t)_fit} |
323 |
C_{2}(t) \approx a e^{-t/\tau_\mathrm{short}} + b e^{-t/\tau_\mathrm{middle}} +\ |
324 |
c |
325 |
e^{-t/\tau_\mathrm{long}} + (1-a-b-c) |
326 |
\end{equation} |
327 |
where $\tau_{short}$ corresponds to the librational motion of the water |
328 |
molecules, $\tau_{middle}$ corresponds to jumps between the breaking and |
329 |
making of hydrogen bonds, and $\tau_{long}$ corresponds to the translational |
330 |
motion of the water molecules. The last term in \eqref{C(t)_fit} accounts |
331 |
for the water molecules trapped in the ice which do not experience any |
332 |
long-time orientational decay. |
333 |
|
334 |
In Figures \ref{fig:PyrOrient} and \ref{fig:SPorient} we see the $z$-coordinate |
335 |
profiles for the three decay constants, $\tau_{short}$ (panel a), |
336 |
$\tau_{middle}$ (panel b), |
337 |
and $\tau_{long}$ (panel c) for the pyramidal and secondary prismatic systems |
338 |
respectively. The control experiments (no shear) are shown in black, and |
339 |
an experiment with an imposed momentum flux is shown in red. The vertical |
340 |
dotted line traversing all three panels denotes the midpoint of the |
341 |
interface as determined by the local tetrahedral order parameter fitting. |
342 |
In the liquid regions of both systems, we see that $\tau_{middle}$ and |
343 |
$\tau_{long}$ have approximately consistent values of $3-6$ ps and $30-40$ ps, |
344 |
resepctively, and increase in value as we approach the interface. Conversely, |
345 |
in panel a, we see that $\tau_{short}$ decreases from the liquid value |
346 |
of $72-76$ fs as we approach the interface. We believe this speed up is due to |
347 |
the constrained motion of librations closer to the interface. Both the |
348 |
approximate values for the decays and trends approaching the interface match |
349 |
those reported previously for the basal and prismatic interfaces. |
350 |
|
351 |
As done previously, we have attempted to quantify the distance, $d_{pyramidal}$ |
352 |
and $d_{secondary prismatic}$, from the |
353 |
interface that the deviations from the bulk liquid values begin. This was done |
354 |
by fitting the orientational decay constant $z$-profiles by |
355 |
\begin{equation}\label{tauFit} |
356 |
\tau(z)\approx\tau_{liquid}+(\tau_{wall}-\tau_{liquid})e^{-(z-z_{wall})/d} |
357 |
\end{equation} |
358 |
where $\tau_{liquid}$ and $\tau_{wall}$ are the liquid and projected wall |
359 |
values of the decay constants, $z_{wall}$ is the location of the interface, |
360 |
and $d$ is the displacement from the interface at which these deviations |
361 |
occur. The values for $d_{pyramidal}$ and $d_{secondary prismatic}$ were |
362 |
determined |
363 |
for each of the decay constants, and then averaged for better statistics |
364 |
($\tau_{middle}$ was ommitted for secondary prismatic). For the pyramidal |
365 |
system, |
366 |
$d_{pyramidal}$ was found to be 2.7 \AA\ for both the control and the sheared |
367 |
system. We found $d_{secondary prismatic}$ to be slightly larger than |
368 |
$d_{pyramidal}$ for both the control and with an applied shear, with |
369 |
displacements of $4$ \AA\ for the control system and $3$ \AA\ for the |
370 |
experiment with the imposed momentum flux. These values are consistent with |
371 |
those found for the basal ($d_{basal}\approx2.9$ \AA\ ) and prismatic |
372 |
($d_{prismatic}\approx3.5$ \AA\ ) systems. |
373 |
|
374 |
\subsection{Coefficient of friction of the interfaces} |
375 |
While investigating the kinetic coefficient of friction, there was found |
376 |
to be a dependence for $\mu_k$ |
377 |
on the temperature of the liquid water in the system. We believe this |
378 |
dependence |
379 |
arrises from the sharp discontinuity of the viscosity for the SPC/E model |
380 |
at temperatures approaching 200 K\cite{kuang12}. Due to this, we propose |
381 |
a weighting to the interfacial friction coefficient, $\kappa$ by the |
382 |
shear viscosity of the fluid at 225 K. The interfacial friction coefficient |
383 |
relates the shear stress with the relative velocity of the fluid normal to the |
384 |
interface: |
385 |
\begin{equation}\label{Shenyu-13} |
386 |
j_{z}(p_{x}) = \kappa[v_{x}(fluid)-v_{x}(solid)] |
387 |
\end{equation} |
388 |
where $j_{z}(p_{x})$ is the applied momentum flux (shear stress) across $z$ |
389 |
in the |
390 |
$x$-dimension, and $v_{x}$(fluid) and $v_{x}$(solid) are the velocities |
391 |
directly adjacent to the interface. The shear viscosity, $\eta(T)$, of the |
392 |
fluid can be determined under a linear response of the momentum |
393 |
gradient to the applied shear stress by |
394 |
\begin{equation}\label{Shenyu-11} |
395 |
j_{z}(p_{x}) = \eta(T) \frac{\partial v_{x}}{\partial z}. |
396 |
\end{equation} |
397 |
Using eqs \eqref{Shenyu-13} and \eqref{Shenyu-11}, we can find the following |
398 |
expression for $\kappa$, |
399 |
\begin{equation}\label{kappa-1} |
400 |
\kappa = \eta(T) \frac{\partial v_{x}}{\partial z}\frac{1}{[v_{x}(fluid)-v_{x}(solid)]}. |
401 |
\end{equation} |
402 |
Here is where we will introduce the weighting term of $\eta(225)/\eta(T)$ |
403 |
giving us |
404 |
\begin{equation}\label{kappa-2} |
405 |
\kappa = \frac{\eta(225)}{[v_{x}(fluid)-v_{x}(solid)]}\frac{\partial v_{x}}{\partial z}. |
406 |
\end{equation} |
407 |
|
408 |
To obtain the value of $\eta(225)$ for the SPC/E model, a $31.09 \times 29.38 |
409 |
\times 124.39$ \AA\ box with 3744 SPC/E liquid water molecules was |
410 |
equilibrated to 225K, |
411 |
and 5 unique shearing experiments were performed. Each experiment was |
412 |
conducted in the NVE and were 5 ns in |
413 |
length. The VSS were attempted every timestep, which was set to 2 fs. |
414 |
For our SPC/E systems, we found $\eta(225)$ to be 0.0148 $\pm$ 0.0007 Pa s, |
415 |
roughly ten times larger than the value found for 280 K SPC/E bulk water by |
416 |
Kuang\cite{kuang12}. |
417 |
|
418 |
The interfacial friction coefficient, $\kappa$, can equivalently be expressed |
419 |
as the ratio of the viscosity of the fluid to the slip length, $\delta$, which |
420 |
is an indication of how 'slippery' the interface is. |
421 |
\begin{equation}\label{kappa-3} |
422 |
\kappa = \frac{\eta}{\delta} |
423 |
\end{equation} |
424 |
In each of the systems, the interfacial temperature was kept fixed to 225K, |
425 |
which ensured the viscosity of the fluid at the |
426 |
interace was approximately the same. Thus, any significant variation in |
427 |
$\kappa$ between |
428 |
the systems indicates differences in the 'slipperiness' of the interfaces. |
429 |
As each of the ice systems are sheared relative to liquid water, the |
430 |
'slipperiness' of the interface can be taken as an indication of how |
431 |
hydrophobic or hydrophilic the interface is. The calculated $\kappa$ values |
432 |
found for the four crystal facets of Ice-I$_\mathrm{h}$ investigated are shown |
433 |
in Table \ref{tab:kapa}. The basal and pyramidal facets were found to have |
434 |
similar values of $\kappa \approx$ 0.0006 |
435 |
(amu \AA\textsuperscript{-2} fs\textsuperscript{-1}), while values of |
436 |
$\kappa \approx$ 0.0003 (amu \AA\textsuperscript{-2} fs\textsuperscript{-1}) |
437 |
were found for the prismatic and secondary prismatic systems. |
438 |
These results indicate that the basal and pyramidal facets are |
439 |
more hydrophilic than the prismatic and secondary prismatic facets. |
440 |
%This indicates something about the similarity between the two facets that |
441 |
%share similar values... |
442 |
%Maybe find values for kappa for other materials to compare against? |
443 |
|
444 |
%\begin{table}[h] |
445 |
%\centering |
446 |
%\caption{$\kappa$ values for the basal, prismatic, pyramidal, and secondary \ |
447 |
%prismatic facets of Ice-I$_\mathrm{h}$} |
448 |
%\label{tab:kappa} |
449 |
%\begin{tabular}{|ccc|} \hline |
450 |
% & \multicolumn{2}{c|}{$\kappa_{Drag direction}$ (amu \AA\textsuperscript{-2} fs\textsuperscript{-1})} \\ |
451 |
% Interface & $\kappa_{x}$ & $\kappa_{y}$ \\ \hline |
452 |
% basal & $0.00059 \pm 0.00003$ & $0.00065 \pm 0.00008$ \\ |
453 |
% prismatic & $0.00030 \pm 0.00002$ & $0.00030 \pm 0.00001$ \\ |
454 |
% pyramidal & $0.00058 \pm 0.00004$ & $0.00061 \pm 0.00005$ \\ |
455 |
% secondary prismatic & $0.00035 \pm 0.00001$ & $0.00033 \pm 0.00002$ \\ \hline |
456 |
%\end{tabular} |
457 |
%\end{table} |
458 |
|
459 |
|
460 |
|
461 |
|
462 |
\begin{table}[h] |
463 |
\centering |
464 |
\caption{Phyiscal properties of the basal, prismatic, pyramidal, and secondary prismatic facets of Ice-I$_\mathrm{h}$} |
465 |
\label{tab:kappa} |
466 |
\begin{tabular}{|ccccc|} \hline |
467 |
& \multicolumn{2}{c}{$\kappa_{Drag direction}$ (amu \AA\textsuperscript{-2} fs\textsuperscript{-1})} & & \\ |
468 |
Interface & $\kappa_{x}$ & $\kappa_{y}$ & $\theta_{\infty}$ & K_{spread} (ns^{-1}) \\ \hline |
469 |
basal & $0.00059 \pm 0.00003$ & $0.00065 \pm 0.00008$ & $34.1 \pm 0.9$ & $0.60 \pm 0.07$ \\ |
470 |
pyramidal & $0.00058 \pm 0.00004$ & $0.00061 \pm 0.00005$ & $35 \pm 3$ & $0.7 \pm 0.1$ \\ |
471 |
prismatic & $0.00030 \pm 0.00002$ & $0.00030 \pm 0.00001$ & $46 \pm 4$ & $0.75 \pm 0.09$ \\ |
472 |
secondary prismatic & $0.00035 \pm 0.00001$ & $0.00033 \pm 0.00002$ & $42 \pm 2$ & $0.69 \pm 0.03$ \\ \hline |
473 |
\end{tabular} |
474 |
\end{table} |
475 |
|
476 |
|
477 |
%\begin{table}[h] |
478 |
%\centering |
479 |
%\caption{Solid-liquid friction coefficients (measured in amu~fs\textsuperscript\ |
480 |
%{-1}). \\ |
481 |
%\textsuperscript{a} See ref. \onlinecite{Louden13}. } |
482 |
%\label{tab:lambda} |
483 |
%\begin{tabular}{|ccc|} \hline |
484 |
% & \multicolumn{2}{c|}{Drag direction} \\ |
485 |
% Interface & $x$ & $y$ \\ \hline |
486 |
% basal\textsuperscript{a} & $0.08 \pm 0.02$ & $0.09 \pm 0.03$ \\ |
487 |
% prismatic (T = 225)\textsuperscript{a} & $0.037 \pm 0.008$ & $0.04 \pm 0.01$ \\ |
488 |
% prismatic (T = 230) & $0.10 \pm 0.01$ & $0.070 \pm 0.006$\\ |
489 |
% pyramidal & $0.13 \pm 0.03$ & $0.14 \pm 0.03$ \\ |
490 |
% secondary prismatic & $0.13 \pm 0.02$ & $0.12 \pm 0.03$ \\ \hline |
491 |
%\end{tabular} |
492 |
%\end{table} |
493 |
|
494 |
|
495 |
\begin{figure} |
496 |
\includegraphics[width=\linewidth]{Pyr-orient} |
497 |
\caption{\label{fig:PyrOrient} The three decay constants of the |
498 |
orientational time correlation function, $C_2(t)$, for water as a function |
499 |
of distance from the center of the ice slab. The vertical dashed line |
500 |
indicates the edge of the pyramidal ice slab determined by the local order |
501 |
tetrahedral parameter. The control (black circles) and sheared (red squares) |
502 |
experiments were fit by a shifted exponential decay (Eq. 9\cite{Louden13}) |
503 |
shown by the black and red lines respectively. The upper two panels show that |
504 |
translational and hydrogen bond making and breaking events slow down |
505 |
through the interface while approaching the ice slab. The bottom most panel |
506 |
shows the librational motion of the water molecules speeding up approaching |
507 |
the ice block due to the confined region of space allowed for the molecules |
508 |
to move in.} |
509 |
\end{figure} |
510 |
|
511 |
\begin{figure} |
512 |
\includegraphics[width=\linewidth]{SP-orient-less} |
513 |
\caption{\label{fig:SPorient} Decay constants for $C_2(t)$ at the secondary |
514 |
prismatic face. Panel descriptions match those in \ref{fig:PyrOrient}.} |
515 |
\end{figure} |
516 |
|
517 |
|
518 |
|
519 |
\section{Conclusion} |
520 |
We present the results of molecular dynamics simulations of the pyrmaidal |
521 |
and secondary prismatic facets of an SPC/E model of the |
522 |
Ice-I$_\mathrm{h}$/water interface. The ice was sheared through the liquid |
523 |
water while being exposed to a thermal gradient to maintain a stable |
524 |
interface by using the minimally perturbing VSS RNEMD method. In agreement |
525 |
with our previous findings for the basal and prismatic facets, the interfacial |
526 |
width was found to be independent of shear rate as measured by the local |
527 |
tetrahedral order parameter. This width was found to be |
528 |
3.2~$\pm$ 0.2~\AA\ for both the pyramidal and the secondary prismatic systems. |
529 |
These values are in good agreement with our previously calculated interfacial |
530 |
widths for the basal (3.2~$\pm$ 0.4~\AA\ ) and prismatic (3.6~$\pm$ 0.2~\AA\ ) |
531 |
systems. |
532 |
|
533 |
Orientational dynamics of the Ice-I$_\mathrm{h}$/water interfaces were studied |
534 |
by calculation of the orientational time correlation function at varying |
535 |
displacements normal to the interface. The decays were fit |
536 |
to a tri-exponential decay, where the three decay constants correspond to |
537 |
the librational motion of the molecules driven by the restoring forces of |
538 |
existing hydrogen bonds ($\tau_{short}$ $\mathcal{O}$(10 fs)), jumps between |
539 |
two different hydrogen bonds ($\tau_{middle}$ $\mathcal{O}$(1 ps)), and |
540 |
translational motion of the molecules ($\tau_{long}$ $\mathcal{O}$(100 ps)). |
541 |
$\tau_{short}$ was found to decrease approaching the interface due to the |
542 |
constrained motion of the molecules as the local environment becomes more |
543 |
ice-like. Conversely, the two longer-time decay constants were found to |
544 |
increase at small displacements from the interface. As seen in our previous |
545 |
work on the basal and prismatic facets, there appears to be a dynamic |
546 |
interface width at which deviations from the bulk liquid values occur. |
547 |
We had previously found $d_{basal}$ and $d_{prismatic}$ to be approximately |
548 |
2.8~\AA\ and 3.5~\AA. We found good agreement of these values for the |
549 |
pyramidal and secondary prismatic systems with $d_{pyramidal}$ and |
550 |
$d_{secondary prismatic}$ to be 2.7~\AA\ and 3~\AA. For all four of the |
551 |
facets, no apparent dependence of the dynamic width on the shear rate was |
552 |
found. |
553 |
|
554 |
%Paragraph summarizing the Kappa values |
555 |
The interfacial friction coefficient, $\kappa$, was determined for each facet |
556 |
interface. We were able to reach an expression for $\kappa$ as a function of |
557 |
the velocity profile of the system which is scaled by the viscosity of the liquid |
558 |
at 225 K. In doing so, we have obtained an expression for $\kappa$ which is |
559 |
independent of temperature differences of the liquid water at far displacements |
560 |
from the interface. We found the basal and pyramidal facets to have |
561 |
similar $\kappa$ values, of $\kappa \approx$ 0.0006 |
562 |
(amu \AA\textsuperscript{-2} fs\textsuperscript{-1}). However, the |
563 |
prismatic and secondary prismatic facets were found to have $\kappa$ values of |
564 |
$\kappa \approx$ 0.0003 (amu \AA\textsuperscript{-2} fs\textsuperscript{-1}). |
565 |
As these ice facets are being sheared relative to liquid water, with the |
566 |
structural and dynamic width of all four |
567 |
interfaces being approximately the same, the difference in the coefficient of |
568 |
friction indicates the hydrophilicity of the crystal facets are not |
569 |
equivalent. Namely, that the basal and pyramidal facets of Ice-I$_\mathrm{h}$ |
570 |
are more hydrophilic than the prismatic and secondary prismatic facets. |
571 |
|
572 |
|
573 |
\begin{acknowledgments} |
574 |
Support for this project was provided by the National |
575 |
Science Foundation under grant CHE-1362211. Computational time was |
576 |
provided by the Center for Research Computing (CRC) at the |
577 |
University of Notre Dame. |
578 |
\end{acknowledgments} |
579 |
|
580 |
\newpage |
581 |
|
582 |
\bibliography{iceWater} |
583 |
|
584 |
\end{document} |