| 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 |
Abstract abstract abstract... |
| 54 |
\end{abstract} |
| 55 |
|
| 56 |
\maketitle |
| 57 |
|
| 58 |
\section{Introduction} |
| 59 |
Explain a little bit about ice Ih, point group stuff. |
| 60 |
|
| 61 |
Mention previous work done / on going work by other people. Haymet and Rick |
| 62 |
seem to be investigating how the interfaces is perturbed by the presence of |
| 63 |
ions. This is the conlcusion of a recent publication of the basal and |
| 64 |
prismatic facets of ice Ih, now presenting the pyramidal and secondary |
| 65 |
prism facets under shear. |
| 66 |
|
| 67 |
\section{Methodology} |
| 68 |
|
| 69 |
\begin{figure} |
| 70 |
\includegraphics[width=\linewidth]{SP_comic_strip} |
| 71 |
\caption{\label{fig:spComic} The secondary prism interface with a shear |
| 72 |
rate of 3.5 ms\textsuperscript{-1}. Lower panel: the local tetrahedral order |
| 73 |
parameter, $q(z)$, (black circles) and the hyperbolic tangent fit (red line). |
| 74 |
Middle panel: the imposed thermal gradient required to maintain a fixed |
| 75 |
interfacial temperature. Upper panel: the transverse velocity gradient that |
| 76 |
develops in response to an imposed momentum flux. The vertical dotted lines |
| 77 |
indicate the locations of the midpoints of the two interfaces.} |
| 78 |
\end{figure} |
| 79 |
|
| 80 |
\begin{figure} |
| 81 |
\includegraphics[width=\linewidth]{Pyr_comic_strip} |
| 82 |
\caption{\label{fig:pyrComic} The pyramidal interface with a shear rate of 3.8 \ |
| 83 |
ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:spComic}.} |
| 84 |
\end{figure} |
| 85 |
|
| 86 |
\subsection{Pyramidal and secondary prism system construction} |
| 87 |
|
| 88 |
The construction of the pyramidal and secondary prism systems follows that of |
| 89 |
the basal and prismatic systems presented elsewhere\cite{Louden13}, however |
| 90 |
the ice crystals and water boxes were equilibrated and combined at 50K |
| 91 |
instead of 225K. The ice / water systems generated were then equilibrated |
| 92 |
to 225K. The resulting pyramidal system was |
| 93 |
$37.47 \times 29.50 \times 93.02$ \AA\ with 1216 |
| 94 |
SPC/E molecules in the ice slab, and 2203 in the liquid phase. The secondary |
| 95 |
prism system generated was $71.87 \times 31.66 \times 161.55$ \AA\ with 3840 |
| 96 |
SPC/E molecules in the ice slab and 8176 molecules in the liquid phase. |
| 97 |
|
| 98 |
\subsection{Computational details} |
| 99 |
% Do we need to justify the sims at 225K? |
| 100 |
% No crystal growth or shrinkage over 2 successive 1 ns NVT simulations for |
| 101 |
% either the pyramidal or sec. prism ice/water systems. |
| 102 |
|
| 103 |
The computational details performed here were equivalent to those reported |
| 104 |
in the previous publication\cite{Louden13}. The only changes made to the |
| 105 |
previously reported procedure were the following. VSS-RNEMD moves were |
| 106 |
attempted every 2 fs instead of every 50 fs. This was done to minimize |
| 107 |
the magnitude of each individual VSS-RNEMD perturbation to the system. |
| 108 |
|
| 109 |
All pyramidal simulations were performed under the NVT ensamble except those |
| 110 |
during which statistics were accumulated for the orientational correlation |
| 111 |
function, which were performed under the NVE ensamble. All secondary prism |
| 112 |
simulations were performed under the NVE ensamble. |
| 113 |
|
| 114 |
\section{Results and discussion} |
| 115 |
\subsection{Interfacial width} |
| 116 |
In the literature there is good agreement that between the solid ice and |
| 117 |
the bulk water, there exists a region of 'slush-like' water molecules. |
| 118 |
In this region, the water molecules are structurely distinguishable and |
| 119 |
behave differently than those of the solid ice or the bulk water. |
| 120 |
The characteristics of this region have been defined by both structural |
| 121 |
and dynamic properties; and its width has been measured by the change of these |
| 122 |
properties from their bulk liquid values to those of the solid ice. |
| 123 |
Examples of these properties include the density, the diffusion constant, and |
| 124 |
the translational order profile. \cite{Bryk02,Karim90,Gay02,Hayward01,Hayward02,Karim88} |
| 125 |
|
| 126 |
Since the VSS-RNEMD moves used to impose the thermal and velocity gradients |
| 127 |
perturb the momenta of the water molecules in |
| 128 |
the systems, parameters that depend on translational motion may give |
| 129 |
faulty results. A stuructural parameter will be less effected by the |
| 130 |
VSS-RNEMD perturbations to the system. Due to this, we have used the |
| 131 |
local order tetrahedral parameter to quantify the width of the interface, |
| 132 |
which was originally described by Kumar\cite{Kumar09} and |
| 133 |
Errington\cite{Errington01} and explained in our |
| 134 |
previous publication\cite{Louden13} in relation to an ice/water system. |
| 135 |
|
| 136 |
Paragraph and eq. for tetrahedrality here. |
| 137 |
|
| 138 |
To determine the width of the interfaces, each of the systems were |
| 139 |
divided into 100 artificial bins along the |
| 140 |
$z$-dimension, and the local tetrahedral order parameter, $q(z)$, was |
| 141 |
time-averaged for each of the bins, resulting in a tetrahedrality profile of |
| 142 |
the system. These profiles are shown across the $z$-dimension of the systems |
| 143 |
in panel $a$ of Figures \ref{fig:spComic} |
| 144 |
and \ref{fig:pyrComic} (black circles). The $q(z)$ function has a range of |
| 145 |
(0,1), where a larger value indicates a more tetrahedral environment. |
| 146 |
The $q(z)$ for the bulk liquid was found to be $\approx $ 0.77, while values of |
| 147 |
$\approx $0.92 were more common for the ice. The tetrahedrality profiles were |
| 148 |
fit using a hyperbolic tangent\cite{Louden13} designed to smoothly fit the |
| 149 |
bulk to ice |
| 150 |
transition, while accounting for the thermal influence on the profile by the |
| 151 |
kinetic energy exchanges of the VSS-RNEMD moves. In panels $b$ and $c$, the |
| 152 |
imposed thermal and velocity gradients can be seen. The verticle dotted |
| 153 |
lines traversing all three panels indicate the midpoints of the interface |
| 154 |
as determined by the hyperbolic tangent fit of the tetrahedrality profiles. |
| 155 |
|
| 156 |
From fitting the tetrahedrality profiles for each of the 0.5 nanosecond |
| 157 |
simulations (panel c of Figures \ref{fig:spComic} and \ref{fig:pyrComic}) |
| 158 |
by Eq. 6\cite{Louden13},we find the interfacial width to be |
| 159 |
$3.2 \pm 0.2$ and $3.2 \pm 0.2$ \AA\ for the control system with no applied |
| 160 |
momentum flux for both the pyramidal and secondary prism systems. |
| 161 |
Over the range of shear rates investigated, |
| 162 |
$0.6 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 5.6 \pm 0.4 \mathrm{ms}^{-1}$ for |
| 163 |
the pyramidal system and $0.9 \pm 0.3 \mathrm{ms}^{-1} \rightarrow 5.4 \pm 0.1 |
| 164 |
\mathrm{ms}^{-1}$ for the secondary prism, we found no significant change in |
| 165 |
the interfacial width. This follows our previous findings of the basal and |
| 166 |
prismatic systems, in which the interfacial width was invarient of the |
| 167 |
shear rate of the ice. The interfacial width of the quiescent basal and |
| 168 |
prismatic systems was found to be $3.2 \pm 0.4$ \AA\ and $3.6 \pm 0.2$ \AA\ |
| 169 |
respectively. Over the range of shear rates investigated, $0.6 \pm 0.3 |
| 170 |
\mathrm{ms}^{-1} \rightarrow 5.3 \pm 0.5 \mathrm{ms}^{-1}$ for the basal |
| 171 |
system and $0.9 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 4.5 \pm 0.1 |
| 172 |
\mathrm{ms}^{-1}$ for the prismatic. |
| 173 |
|
| 174 |
These results indicate that the surface structure of the exposed ice crystal |
| 175 |
has little to no effect on how far into the bulk the ice-like structural |
| 176 |
ordering is. Also, it appears that the interface is not structurally effected |
| 177 |
by shearing the ice through water. |
| 178 |
|
| 179 |
|
| 180 |
\subsection{Orientational dynamics} |
| 181 |
%Should we include the math here? |
| 182 |
The orientational time correlation function, |
| 183 |
\begin{equation}\label{C(t)1} |
| 184 |
C_{2}(t)=\langle P_{2}(\mathbf{u}(0)\cdot \mathbf{u}(t))\rangle, |
| 185 |
\end{equation} |
| 186 |
helps indicate the local environment around the water molecules. The function |
| 187 |
begins with an initial value of unity, and decays to zero as the water molecule |
| 188 |
loses memory of its former orientation. Observing the rate at which this decay |
| 189 |
occurs can provide insight to the mechanism and timescales for the relaxation. |
| 190 |
In eq. \eqref{C(t)1}, $P_{2}$ is the second-order Legendre polynomial, and |
| 191 |
$\mathbf{u}$ is the the bisecting HOH vector. The angle brackets indicate |
| 192 |
an ensemble average over all the water molecules in a given spatial region. |
| 193 |
|
| 194 |
To investigate the dynamics of the water molecules across the interface, the |
| 195 |
systems were divided in the $z$-dimension into bins, each $\approx$ 3 \AA\ |
| 196 |
wide, and \eqref{C(t)1} was computed for each of the bins. A water |
| 197 |
molecule was allocated to a particular bin if it was initially in the bin |
| 198 |
at time zero. To compute \eqref{C(t)1}, each 0.5 ns simulation was followed |
| 199 |
by an additional 200 ps microcanonical (NVE) simulation during which the |
| 200 |
position and orientations of each molecule were recorded every 0.1 ps. |
| 201 |
|
| 202 |
The data obtained for each bin was then fit to a triexponential decay given by |
| 203 |
\begin{equation}\label{C(t)_fit} |
| 204 |
C_{2}(t) \approx a e^{-t/\tau_\mathrm{short}} + b e^{-t/\tau_\mathrm{middle}} +\ |
| 205 |
c |
| 206 |
e^{-t/\tau_\mathrm{long}} + (1-a-b-c) |
| 207 |
\end{equation} |
| 208 |
where $\tau_{short}$ corresponds to the librational motion of the water |
| 209 |
molecules, $\tau_{middle}$ corresponds to jumps between the breaking and |
| 210 |
making of hydrogen bonds, and $\tau_{long}$ corresponds to the translational |
| 211 |
motion of the water molecules. The last term in \eqref{C(t)_fit} accounts |
| 212 |
for the water molecules trapped in the ice which do not experience any |
| 213 |
long-time orientational decay. |
| 214 |
|
| 215 |
In Figures \ref{fig:PyrOrient} and \ref{fig:SPorient} we see the $z$-coordinate |
| 216 |
profiles for the three decay constants, $\tau_{short}$ (panel a), |
| 217 |
$\tau_{middle}$ (panel b), |
| 218 |
and $\tau_{long}$ (panel c) for the pyramidal and secondary prismatic systems |
| 219 |
respectively. The control experiments (no shear) are shown in black, and |
| 220 |
an experiment with an imposed momentum flux is shown in red. The vertical |
| 221 |
dotted line traversing all three panels denotes the midpoint of the |
| 222 |
interface as determined by the local tetrahedral order parameter fitting. |
| 223 |
In the liquid regions of both systems, we see that $\tau_{middle}$ and |
| 224 |
$\tau_{long}$ have approximately consistent values of $3-6$ ps and $30-40$ ps, |
| 225 |
resepctively, and increase in value as we approach the interface. Conversely, |
| 226 |
in panel a, we see that $\tau_{short}$ decreases from the liquid value |
| 227 |
of $72-76$ fs as we approach the interface. We believe this speed up is due to |
| 228 |
the constrained motion of librations closer to the interface. Both the |
| 229 |
approximate values for the decays and relative trends match those reported |
| 230 |
previously for the basal and prismatic interfaces. |
| 231 |
|
| 232 |
As done previously, we have attempted to quantify the distance, $d_{pyramidal}$ |
| 233 |
and $d_{secondary prism}$, from the |
| 234 |
interface that the deviations from the bulk liquid values begin. This was done |
| 235 |
by fitting the orientational decay constant $z$-profiles by |
| 236 |
\begin{equation}\label{tauFit} |
| 237 |
\tau(z)\approx\tau_{liquid}+(\tau_{solid}-\tau_{liquid})e^{-(z-z_{wall})/d} |
| 238 |
\end{equation} |
| 239 |
where $\tau_{liquid}$ and $\tau_{solid}$ are the liquid and projected solid |
| 240 |
values of the decay constants, $z_{wall}$ is the location of the interface, |
| 241 |
and $d$ is the displacement from the interface at which these deviations |
| 242 |
occur. The values for $d_{pyramidal}$ and $d_{secondary prismatic}$ were |
| 243 |
determined |
| 244 |
for each of the decay constants, and then averaged for better statistics |
| 245 |
($\tau_{middle}$ was ommitted for secondary prism). For the pyramidal system, |
| 246 |
$d_{pyramidal}$ was found to be 2.7 \AA\ for both the control and the sheared |
| 247 |
system. We found $d_{secondary prismatic}$ to be slightly larger than |
| 248 |
$d_{pyramidal}$ for both the control and with an applied shear, with |
| 249 |
displacements of $4$ \AA\ for the control system and $3$ \AA\ for the |
| 250 |
experiment with the imposed momentum flux. These values are consistent with |
| 251 |
those found for the basal ($d_{basal}\approx2.9$ \AA\ ) and prismatic |
| 252 |
($d_{prismatic}\approx3.5$ \AA\ ) systems. |
| 253 |
|
| 254 |
\subsection{Coefficient of friction of the interfaces} |
| 255 |
While investigating the kinetic coefficient of friction for the larger |
| 256 |
prismatic system, there was found to be a dependence for $\mu_k$ |
| 257 |
on the temperature of the liquid water in the system. We believe this |
| 258 |
dependence |
| 259 |
arrises from the sharp discontinuity of the viscosity for the SPC/E model |
| 260 |
at temperatures approaching 200 K\cite{kuang12}. Due to this, we propose |
| 261 |
a weighting to the structural interfacial parameter, $\kappa$ by the |
| 262 |
viscosity at $225$ K, the temperature of the interface. $\kappa$ is |
| 263 |
traditionally defined as |
| 264 |
\begin{equation}\label{kappa} |
| 265 |
\kappa = \eta/\delta |
| 266 |
\end{equation} |
| 267 |
where $\eta$ is the viscosity and $\delta$ is the slip length. |
| 268 |
In our ice/water shearing simulations, the system has reached a steady state |
| 269 |
when the applied force, |
| 270 |
|
| 271 |
\begin{equation} |
| 272 |
f_{applied} = \mathbf{j}_z(\mathbf{p})L_x L_y |
| 273 |
\end{equation} |
| 274 |
is equal to the frictional force resisting the motion of the ice block |
| 275 |
\begin{equation} |
| 276 |
f_{friction} = \frac{\eta \mathbf{v} L_x L_y}{\delta} |
| 277 |
\end{equation} |
| 278 |
where $\mathbf{v}$ is the relative velocity of the liquid from the ice. |
| 279 |
When this condition is met, we are able to solve the resulting expression to |
| 280 |
obtain, |
| 281 |
\begin{equation}\label{force_equality} |
| 282 |
\frac{\eta}{\delta} = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}} |
| 283 |
\end{equation} |
| 284 |
From \eqref{kappa}, \eqref{force_equality} becomes |
| 285 |
\begin{equation} |
| 286 |
\kappa = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}} |
| 287 |
\end{equation} |
| 288 |
which we will multiply by a viscosity weighting term to reach |
| 289 |
\begin{equation} \label{kappa2} |
| 290 |
\kappa = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}} \frac{\eta(225)}{\eta(T)} |
| 291 |
\end{equation} |
| 292 |
Assuming linear response theory is valid, an expression for ($\eta$) can |
| 293 |
be found from the imposed momentum flux and the measured velocity gradient. |
| 294 |
\begin{equation}\label{eta_eq} |
| 295 |
\eta = \frac{\mathbf{j}_z(\mathbf{p})}{\frac{\partial v_x}{\partial z}} |
| 296 |
\end{equation} |
| 297 |
Substituting eq \eqref{eta_eq} into eq \eqref{kappa2} we arrive at |
| 298 |
\begin{equation} |
| 299 |
\kappa = \frac{\frac{\partial v_x}{\partial z}}{\mathbf{v}}\eta(225) |
| 300 |
\end{equation} |
| 301 |
|
| 302 |
\begin{table}[h] |
| 303 |
\centering |
| 304 |
\caption{$\kappa$ values for the basal, prismatic, pyramidal, and secondary prismatic facets of Ice-I$_\mathrm{h}$} |
| 305 |
\label{tab:kappa} |
| 306 |
\begin{tabular}{|ccc|} \hline |
| 307 |
& \multicolumn{2}{c|}{$\kappa_{Drag direction}$} \\ |
| 308 |
Interface & $\kappa_{x}$ & $\kappa_{y}$ \\ \hline |
| 309 |
basal & $0.00059 \pm 0.00003$ & $0.00065 \pm 0.00008$ \\ |
| 310 |
prismatic & $0.00030 \pm 0.00002$ & $0.00030 \pm 0.00001$ \\ |
| 311 |
pyramidal & $0.00058 \pm 0.00004$ & $0.00061 \pm 0.00005$ \\ |
| 312 |
secondary prism & $0.00035 \pm 0.00001$ & $0.00033 \pm 0.00002$ \\ \hline |
| 313 |
\end{tabular} |
| 314 |
\end{table} |
| 315 |
|
| 316 |
|
| 317 |
To obtain the value of $\eta(225)$ for the SPC/E model, a $31.09 \times 29.38 |
| 318 |
\times 124.39$ \AA\ box with 3744 water molecules was equilibrated to 225K, |
| 319 |
and 5 unique shearing experiments were performed. Each experiment was |
| 320 |
conducted in the microcanonical ensemble (NVE) and were 5 ns in |
| 321 |
length. The VSS were attempted every timestep, which was set to 2 fs. |
| 322 |
For our SPC/E systems, we found $\eta(225)$ to be $0.0148 \pm 0.0007 Pa s$, |
| 323 |
roughly ten times larger than the value found for 280 K SPC/E water by |
| 324 |
Kuang\cite{kuang12}. |
| 325 |
|
| 326 |
|
| 327 |
\begin{table}[h] |
| 328 |
\centering |
| 329 |
\caption{Solid-liquid friction coefficients (measured in amu~fs\textsuperscript\ |
| 330 |
{-1}). \\ |
| 331 |
\textsuperscript{a} See ref. \onlinecite{Louden13}. } |
| 332 |
\label{tab:lambda} |
| 333 |
\begin{tabular}{|ccc|} \hline |
| 334 |
& \multicolumn{2}{c|}{Drag direction} \\ |
| 335 |
Interface & $x$ & $y$ \\ \hline |
| 336 |
basal\textsuperscript{a} & $0.08 \pm 0.02$ & $0.09 \pm 0.03$ \\ |
| 337 |
prismatic (T = 225)\textsuperscript{a} & $0.037 \pm 0.008$ & $0.04 \pm 0.01$ \\ |
| 338 |
prismatic (T = 230) & $0.10 \pm 0.01$ & $0.070 \pm 0.006$\\ |
| 339 |
pyramidal & $0.13 \pm 0.03$ & $0.14 \pm 0.03$ \\ |
| 340 |
secondary prism & $0.13 \pm 0.02$ & $0.12 \pm 0.03$ \\ \hline |
| 341 |
\end{tabular} |
| 342 |
\end{table} |
| 343 |
|
| 344 |
|
| 345 |
\begin{figure} |
| 346 |
\includegraphics[width=\linewidth]{Pyr-orient} |
| 347 |
\caption{\label{fig:PyrOrient} The three decay constants of the |
| 348 |
orientational time correlation function, $C_2(t)$, for water as a function |
| 349 |
of distance from the center of the ice slab. The vertical dashed line |
| 350 |
indicates the edge of the pyramidal ice slab determined by the local order |
| 351 |
tetrahedral parameter. The control (black circles) and sheared (red squares) |
| 352 |
experiments were fit by a shifted exponential decay (Eq. 9\cite{Louden13}) |
| 353 |
shown by the black and red lines respectively. The upper two panels show that |
| 354 |
translational and hydrogen bond making and breaking events slow down |
| 355 |
through the interface while approaching the ice slab. The bottom most panel |
| 356 |
shows the librational motion of the water molecules speeding up approaching |
| 357 |
the ice block due to the confined region of space allowed for the molecules |
| 358 |
to move in.} |
| 359 |
\end{figure} |
| 360 |
|
| 361 |
\begin{figure} |
| 362 |
\includegraphics[width=\linewidth]{SP-orient-less} |
| 363 |
\caption{\label{fig:SPorient} Decay constants for $C_2(t)$ at the secondary |
| 364 |
prism face. Panel descriptions match those in \ref{fig:PyrOrient}.} |
| 365 |
\end{figure} |
| 366 |
|
| 367 |
|
| 368 |
|
| 369 |
\section{Conclusion} |
| 370 |
Conclude conclude conclude... |
| 371 |
|
| 372 |
|
| 373 |
\begin{acknowledgments} |
| 374 |
Support for this project was provided by the National |
| 375 |
Science Foundation under grant CHE-1362211. Computational time was |
| 376 |
provided by the Center for Research Computing (CRC) at the |
| 377 |
University of Notre Dame. |
| 378 |
\end{acknowledgments} |
| 379 |
|
| 380 |
\newpage |
| 381 |
|
| 382 |
\bibliography{iceWater} |
| 383 |
|
| 384 |
\end{document} |
