Molecular interpretation of the non-Newtonian viscoelastic behavior of liquid water at high frequencies

Using classical as well as ab initio molecular dynamics simulations, we calculate the frequency-dependent shear viscosity of pure water and water–glycerol mixtures. In agreement with recent experiments, we find deviations from Newtonian-fluid behavior in the THz regime. Based on an extension of the Maxwell model, we introduce a viscoelastic model to describe the observed viscosity spectrum of pure water. We find four relaxation modes in the spectrum which we attribute to (i) hydrogen–bond network topology changes, (ii) hydrogen–bond stretch vibrations of water pairs, (iii) collective vibrations of water molecule triplets, and (iv) librational excitations of individual water molecules. Our model quantitatively describes the viscoelastic response of liquid water on short timescales, where the hydrodynamic description via a Newtonian-fluid model breaks down.

Figure 1

Viscosity spectra $\tilde{\eta }$ as calculated from MD simulations of TIP4P/2005 water at different glycerol concentrations. The real and imaginary parts of the spectra (blue and green solid lines) are calculated from the MD data at $\mathbf{k}=0$ using Eq. (8). The result of a Maxwell model fit, cf., Eq. (13), to frequencies ${10}^{-4}\mathrm{THz}<f<1\mathrm{THz}$, is shown as orange and violet dashed lines.

Figure 2

Comparison of experimental and simulation data for the viscoelasticity of glycerol solutions. MD data points are obtained using the Green-Kubo formula Eq. (9) (red crosses), or fitting a Maxwell model Eq. (13) to the full viscosity spectra (blue dots), cf., Fig. 1. Experimental data for ${\eta }_0$ are taken from Cheng [30], timescales are calculated using $\tau ={\eta }_0/K$, with ${\eta }_0$ taken from Cheng [30] and $K$ taken from Slie et al. [6]

Figure 3

Mean-squared displacement (MSD) of individual water molecules in glycerol solutions. (a) For each glycerol mass fraction, the MSD of each water molecule center of mass is calculated and an average over all such MSDs is performed (colored solid lines). The vertical dashed lines denote the corresponding fitted Maxwell-model timescales $\tau$ from Fig. 2, the horizontal gray shaded region denotes the range $[0.2^2{\mathrm{nm}}^2,0.5^2{\mathrm{nm}}^2]$. The long-time diffusion coefficient $D_{\infty }$, defined by $\mathrm{MSD}=6D_{\infty }t$ for $t$ large enough, is obtained from each curve by fitting a linear function to each dataset for the time interval $t/\mathrm{ps}\in [5\times {10}^4,{10}^5]$. (b) The colored lines show the MSDs from subplot (a) divided by the respective long-time diffusive MSD, and with time rescaled by the respective Maxwell timescale $\tau$. The time $t=\tau$ is indicated by a vertical solid line, the curve $\mathrm{MSD}=6D_{\infty }t$ is depicted as a horizontal solid line.

Figure 4

Complex viscosities of viscoelastic materials. (a) Complex viscosity of a Kelvin–Voigt model, $\tilde{\eta }\left(\omega \right)={\eta }_0+K/\left(-i\omega \right)$, with ${\eta }_0=1\mathrm{mPa}\mathrm{s}$, $K=5\times {10}^{10}\mathrm{mPa}$. (b) Complex viscosity of a Maxwell model, Eq. (13), with ${\eta }_0=1\mathrm{mPa}\mathrm{s}$, $K=2\times {10}^{12}\mathrm{mPa}$. (c) Complex viscosity of a shear inertia model, Eq. (16), with ${\eta }_0=1\mathrm{mPa}\mathrm{s}$, $K=2\times {10}^{12}\mathrm{mPa}\mathrm{s}$, $L={10}^{-12}\mathrm{mPa}{\mathrm{s}}^2$.

Figure 5

Viscosity spectrum of pure water. (a): Complex viscosity $\tilde{\eta }$ from MD simulations of TIP4P/2005 water, calculated using Eq. (8), together with a fit of the viscoelastic network shown in subplot (d) with complex viscosity Eq. (18). The fitting parameters are given in Table 3. (b), (e): Real and imaginary parts of the individual constituents of the fitted viscosity from subplot (a); (c): Complex viscosity $\tilde{\eta }$ from MD simulations of TIP4P/2005 water (replot of subplot (a)), together with the viscosity calculated from ab initio molecular dynamics (aiMD) simulations using Eq. (10) at finite $k=4.016\mathrm{nm}{}^{-1}$. (d): Viscoelastic network used for modeling the shear response of pure water. The viscoelastic network is a parallel combination of two shear inertia models and two Maxwell models, the resulting complex viscosity is given by Eq. (18); (f): Black line: Potential of mean force $U\left(r\right)$ for the oxygen-oxygen distance in TIP4P/2005 water. Blue and yellow lines: Harmonic potentials, fitted around the first and second minima of the pmf. The colored arrows indicate the microscopic interpretation of the constituents of the viscoelastic model shown in subplot (d), namely escape from the nearest neighbor pmf minimum (I), oscillations around the minima (II, III) and librations of individual water molecules (IV).