| 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 24 \AA\ x 36 \AA\ x 99 \AA\ with 900 SPC/E molecules in the ice slab and 1846 SPC/E molecules in the liquid phase. Similarly, the prismatic system was constructed as 36 \AA\ x 36 \AA\ x 86 \AA\ with 1000 SPC/E molecules in the ice slab and 2684 SPC/E molecules in the liquid phase. |
| 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. |
| 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. |
| 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 |
|
|
| 366 |
|
\section{Results and discussion} |
| 367 |
|
|
| 379 |
|
influenced by the RNEMD perturbations to the same degree. Here, we |
| 380 |
|
have used the local tetrahedral order parameter as described by |
| 381 |
|
Kumar\cite{Kumar09} and Errington\cite{Errington01} as our principal |
| 382 |
< |
measure of the interfacial width. |
| 382 |
> |
measure of the interfacial width. A previous study by Bryk and Haymet |
| 383 |
> |
also used local tetrahedrality as an order parameter for ice/water |
| 384 |
> |
interfaces.\cite{Bryk2004b} |
| 385 |
|
|
| 386 |
|
The local tetrahedral order parameter, $q(z)$, is given by |
| 387 |
|
\begin{equation} |
| 393 |
|
where $\psi_{ikj}$ is the angle formed between the oxygen site on |
| 394 |
|
central molecule $k$, and the oxygen sites on two of the four closest |
| 395 |
|
molecules, $i$ and $j$. Molecules $i$ and $j$ are further restricted |
| 396 |
< |
to lie within the first solvation shell of molecule $k$. $N_z = \int |
| 396 |
> |
to lie withing the first peak in the pair distribution function for |
| 397 |
> |
molecule $k$ (typically $<$ 3.41 \AA\ for water). $N_z = \int |
| 398 |
|
\delta(z_k - z) \mathrm{d}z$ is a normalization factor to account for |
| 399 |
|
the varying population of molecules within each finite-width bin. The |
| 400 |
|
local tetrahedral order parameter has a range of $(0,1)$, where the |
| 405 |
|
$q(z) \approx 0.75$ are more common. |
| 406 |
|
|
| 407 |
|
To estimate the interfacial width, the system was divided into 100 |
| 408 |
< |
bins along the $z$-dimension, and a cutoff radius for the first |
| 409 |
< |
solvation shell was set to 3.41~\AA\ . The $q_{z}$ function was |
| 410 |
< |
time-averaged to give yield a tetrahedrality profile of the |
| 411 |
< |
system. The profile was then fit to a hyperbolic tangent that smoothly |
| 403 |
< |
links the liquid and solid states, |
| 408 |
> |
bins along the $z$-dimension. The $q_{z}$ function was time-averaged |
| 409 |
> |
to give yield a tetrahedrality profile of the system. The profile was |
| 410 |
> |
then fit to a hyperbolic tangent that smoothly links the liquid and |
| 411 |
> |
solid states, |
| 412 |
|
\begin{equation}\label{tet_fit} |
| 413 |
|
q(z) \approx |
| 414 |
|
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- |
| 419 |
|
the interface. $l$ and $r$ are the midpoints of the left and right |
| 420 |
|
interfaces, respectively. The last term in eq. \eqref{tet_fit} |
| 421 |
|
accounts for the influence that the weak thermal gradient has on the |
| 422 |
< |
tetrahedrality profile in the liquid region. To estimate the |
| 415 |
< |
10\%-90\% widths commonly used in previous studies,\cite{Bryk02} it is |
| 416 |
< |
a simple matter to scale the widths obtained from the hyperbolic |
| 417 |
< |
tangent fits to obtain $w_{10-90} = 2.1971 \times w$.\cite{Bryk02} |
| 422 |
> |
tetrahedrality profile in the liquid region. |
| 423 |
|
|
| 424 |
|
In Figures \ref{fig:bComic} and \ref{fig:pComic} we see the |
| 425 |
|
$z$-coordinate profiles for tetrahedrality, temperature, and the |
| 434 |
|
225~$\pm$~5K, can be seen in middle panels. The transverse velocity |
| 435 |
|
profile is shown in the upper panels. It is clear from the upper |
| 436 |
|
panels that water molecules in close proximity to the surface (i.e. |
| 437 |
< |
within 10~\AA\ to 15~\AA~of the interfaces) have transverse |
| 438 |
< |
velocities quite close to the velocities within the ice block. There |
| 439 |
< |
is no velocity discontinuity at the interface, which indicates that |
| 440 |
< |
the shearing of ice/water interfaces occurs in the ``stick'' or |
| 441 |
< |
no-slip boundary conditions. |
| 437 |
> |
within 10~\AA\ to 15~\AA~of the interfaces) have transverse velocities |
| 438 |
> |
quite close to the velocities within the ice block. There is no |
| 439 |
> |
velocity discontinuity at the interface, which indicates that the |
| 440 |
> |
shearing of ice/water interfaces occurs in the ``stick'' or no-slip |
| 441 |
> |
boundary conditions. |
| 442 |
|
|
| 443 |
|
\begin{figure} |
| 444 |
|
\includegraphics[width=\linewidth]{bComicStrip} |
| 445 |
< |
\caption{\label{fig:bComic} The basal interfaces. Lower panel: the |
| 446 |
< |
local tetrahedral order parameter, $q(z)$, (black circles) and the |
| 447 |
< |
hyperbolic tangent fit (red line). Middle panel: the imposed |
| 448 |
< |
thermal gradient required to maintain a fixed interfacial |
| 449 |
< |
temperature. Upper panel: the transverse velocity gradient that |
| 450 |
< |
develops in response to an imposed momentum flux. The vertical |
| 451 |
< |
dotted lines indicate the locations of the midpoints of the two |
| 452 |
< |
interfaces.} |
| 445 |
> |
\caption{\label{fig:bComic} The basal interface with a shear rate of |
| 446 |
> |
XXXX. Lower panel: the local tetrahedral order parameter, $q(z)$, |
| 447 |
> |
(black circles) and the hyperbolic tangent fit (red line). Middle |
| 448 |
> |
panel: the imposed thermal gradient required to maintain a fixed |
| 449 |
> |
interfacial temperature. Upper panel: the transverse velocity |
| 450 |
> |
gradient that develops in response to an imposed momentum flux. The |
| 451 |
> |
vertical dotted lines indicate the locations of the midpoints of the |
| 452 |
> |
two interfaces.} |
| 453 |
|
\end{figure} |
| 454 |
|
|
| 455 |
|
\begin{figure} |
| 456 |
|
\includegraphics[width=\linewidth]{pComicStrip} |
| 457 |
< |
\caption{\label{fig:pComic} The prismatic interfaces. Panel |
| 457 |
> |
\caption{\label{fig:pComic} The prismatic interface with a shear rate |
| 458 |
> |
of XXXX. Panel |
| 459 |
|
descriptions match those in figure \ref{fig:bComic}} |
| 460 |
|
\end{figure} |
| 461 |
|
|
| 462 |
< |
From the fits using eq. \eqref{tet_fit}, we find the interfacial |
| 463 |
< |
width for the basal and prismatic systems to be 3.2~$\pm$~0.4~\AA\ and |
| 462 |
> |
From the fits using eq. \eqref{tet_fit}, we find the interfacial width |
| 463 |
> |
for the basal and prismatic systems to be 3.2~$\pm$~0.4~\AA\ and |
| 464 |
|
3.6~$\pm$~0.2~\AA\ , respectively, with no applied momentum flux. Over |
| 465 |
|
the range of shear rates investigated, $0.6 \pm 0.3 \mathrm{ms}^{-1} |
| 466 |
|
\rightarrow 5.3 \pm 0.5 \mathrm{ms}^{-1}$ for the basal system and |
| 468 |
|
\mathrm{ms}^{-1}$ for the prismatic, we found no appreciable change in |
| 469 |
|
the interface width. The fit values for the interfacial width ($w$) |
| 470 |
|
over all shear rates contained the values reported above within their |
| 471 |
< |
error bars. |
| 471 |
> |
error bars. Note that the interfacial widths reported here are based |
| 472 |
> |
on the hyperbolic tangent parameter $w$ in Eq. \ref{tet_fit}. This is |
| 473 |
> |
related to, but not identical with, the 10\%-90\% intefacial widths |
| 474 |
> |
commonly used in previous studies.\cite{Bryk02,Bryk2004b} To estimate |
| 475 |
> |
the 10\%-90\% widths, it is a simple matter to scale the widths |
| 476 |
> |
obtained from the hyperbolic tangent fits to obtain $w_{10-90} = |
| 477 |
> |
2.1971 \times w$.\cite{Bryk02,Bryk2004b} This results in $w_{10-90}$ |
| 478 |
> |
values of 7.0~$\pm$~0.9~\AA\ for the basal face, and 7.9~$\pm$~0.4 |
| 479 |
> |
\AA\ for the prismatic face. These are somewhat smaller than |
| 480 |
> |
previously reported values. |
| 481 |
> |
|
| 482 |
> |
\begin{figure} |
| 483 |
> |
\includegraphics[width=\linewidth]{interface_width_by_shear_rate} |
| 484 |
> |
\caption{\label{fig:widthByShear} The width of the ice water |
| 485 |
> |
interfaces (as measured by Eq. \ref{tet_fit}) exhibits no dependence |
| 486 |
> |
on the applied shear rate between the ice and water regions.} |
| 487 |
> |
\end{figure} |
| 488 |
> |
|
| 489 |
> |
|
| 490 |
|
|
| 491 |
|
\subsubsection{Orientational Dynamics} |
| 492 |
|
The orientational time correlation function, |
| 641 |
|
\lambda=\frac{\eta}{\delta} |
| 642 |
|
\end{equation} |
| 643 |
|
where $\delta$ is the slip length for the liquid measured at the |
| 644 |
< |
location of the interface itself. |
| 644 |
> |
location of the interface itself. In our simulations, the shoulder on |
| 645 |
> |
the velocity profile indicating the location of the hydrodynamic |
| 646 |
> |
boundary in the liquid is not always apparent. In some cases, the |
| 647 |
> |
linear behavior persists nearly up to the interfacial region. For |
| 648 |
> |
this reason, the hydrodynamic position of the boundary is not always |
| 649 |
> |
computable, while the Pit approach (Eq. \ref{Pit}) can be used to find |
| 650 |
> |
the solid-liquid friction coefficient more reliably. |
| 651 |
|
|
| 652 |
< |
In both of these expressions, $\eta$ is the shear viscosity of the |
| 653 |
< |
bulk-like region of the liquid, a quantity which is easily obtained in |
| 654 |
< |
VSS-RNEMD simulations by fitting the velocity profile in the region |
| 655 |
< |
far from the surface.\cite{Kuang12} Assuming linear response in the |
| 656 |
< |
bulk-like region, |
| 652 |
> |
In both the Pit and hydrodynamic boundary expressions, $\eta$ is the |
| 653 |
> |
shear viscosity of the bulk-like region of the liquid, a quantity |
| 654 |
> |
which is easily obtained in VSS-RNEMD simulations by fitting the |
| 655 |
> |
velocity profile in the region far from the surface.\cite{Kuang12} |
| 656 |
> |
Assuming linear response in the bulk-like region, |
| 657 |
|
\begin{equation}\label{Kuang} |
| 658 |
|
j_{z}(p_{x})=-\eta \left(\frac{\partial v_{x}}{\partial z}\right) |
| 659 |
|
\end{equation} |
| 664 |
|
z}\right) \delta} |
| 665 |
|
\end{equation} |
| 666 |
|
|
| 667 |
< |
For ice / water interfaces, the boundary conditions are markedly |
| 668 |
< |
no-slip, so projecting the bulk liquid state velocity profile yields a |
| 669 |
< |
negative slip length. This length is the difference between the |
| 670 |
< |
structural edge of the ice (determined by the tetrahedrality profile) |
| 671 |
< |
and the location where the projected velocity of the bulk liquid |
| 672 |
< |
intersects the solid phase velocity (see Figure |
| 673 |
< |
\ref{fig:delta_example}). The coefficients of friction for the basal |
| 674 |
< |
and the prismatic facets are found to be |
| 675 |
< |
0.07~$\pm$~0.01~amu~fs\textsuperscript{-1} and |
| 646 |
< |
0.032~$\pm$~0.007~amu~fs\textsuperscript{-1}, respectively. These |
| 647 |
< |
results may seem surprising as the basal face is smoother than the |
| 648 |
< |
prismatic with only small undulations of the oxygen positions, while |
| 649 |
< |
the prismatic surface has deep corrugated channels. The applied |
| 650 |
< |
momentum flux used in our simulations is parallel to these channels, |
| 651 |
< |
however, and this results in a flow of water in the same direction as |
| 652 |
< |
the corrugations, allowing water molecules to pass through the |
| 653 |
< |
channels during the shear. |
| 667 |
> |
For ice / water interfaces, the boundary conditions are no-slip, so |
| 668 |
> |
projecting the bulk liquid state velocity profile yields a negative |
| 669 |
> |
slip length. This length is the difference between the structural edge |
| 670 |
> |
of the ice (determined by the tetrahedrality profile) and the location |
| 671 |
> |
where the projected velocity of the bulk liquid intersects the solid |
| 672 |
> |
phase velocity (see Figure \ref{fig:delta_example}). The coefficients |
| 673 |
> |
of friction for the basal and the prismatic facets were determined for |
| 674 |
> |
shearing along both the $x$ and $y$ axes. The values are given in |
| 675 |
> |
table \ref{tab:lambda}. |
| 676 |
|
|
| 677 |
< |
\begin{figure} |
| 677 |
> |
Note that the measured friction coefficient for the basal face is |
| 678 |
> |
twice that of the prismatic face (regardless of drag direction). |
| 679 |
> |
These results may seem surprising as the basalface appears smoother |
| 680 |
> |
than the prismatic with only small undulations of the oxygen |
| 681 |
> |
positions, while the prismatic surface has deep corrugated channels |
| 682 |
> |
along the $x$ direction in the crystal system used in this work. |
| 683 |
> |
However, the corrugations are relatively thin, and the liquid phase |
| 684 |
> |
water does not appear to populate the channels. The prismatic face |
| 685 |
> |
therefore effectively presents stripes of solid-phase molecules |
| 686 |
> |
(making up approximately half of the exposed surface area) with nearly |
| 687 |
> |
empty space between them. The interfacial friction appears to be |
| 688 |
> |
independent of the drag direction, so flow parallel to these channels |
| 689 |
> |
does not explain the lower friction of the prismatic face. A more |
| 690 |
> |
likely explanation is that the effective contact between the liquid |
| 691 |
> |
phase and the prismatic face is reduced by the empty corrugations. |
| 692 |
> |
|
| 693 |
> |
\begin{figure} |
| 694 |
|
\includegraphics[width=\linewidth]{delta_example} |
| 695 |
|
\caption{\label{fig:delta_example} Determining the (negative) slip |
| 696 |
|
length ($\delta$) for the ice-water interfaces (which have decidedly |
| 699 |
|
profile) and the location where the projected velocity of the bulk |
| 700 |
|
liquid (dashed red line) intersects the solid phase velocity (solid |
| 701 |
|
black line). The dotted line indicates the location of the ice as |
| 702 |
< |
determined by the tetrahedrality profile.} |
| 702 |
> |
determined by the tetrahedrality profile. This example is taken |
| 703 |
> |
from the basal-face simulation with an applied shear rate of XXXX.} |
| 704 |
|
\end{figure} |
| 705 |
|
|
| 706 |
|
|
| 707 |
+ |
\begin{table}[h] |
| 708 |
+ |
\centering |
| 709 |
+ |
\caption{Solid-liquid friction coefficients (measured in amu~fs\textsuperscript{-1}) } |
| 710 |
+ |
\label{tab:lambda} |
| 711 |
+ |
\begin{tabular}{|ccc|} \hline |
| 712 |
+ |
& \multicolumn{2}{c|}{Drag direction} \\ |
| 713 |
+ |
Interface & $x$ & $y$ \\ \hline |
| 714 |
+ |
basal & $0.08 \pm 0.02$ & $0.09 \pm 0.03$ \\ |
| 715 |
+ |
prismatic & $0.037 \pm 0.008$ & $0.04 \pm 0.01$ \\ \hline |
| 716 |
+ |
\end{tabular} |
| 717 |
+ |
\end{table} |
| 718 |
+ |
|
| 719 |
+ |
|
| 720 |
|
\section{Conclusion} |
| 721 |
|
We have simulated the basal and prismatic facets of an SPC/E model of |
| 722 |
|
the ice I$_\mathrm{h}$ / water interface. Using VSS-RNEMD, the ice |
