Non-Newtonian flow effects in supercooled water

The viscosity of supercooled water has been a subject of intense study, in particular with respect to its temperature dependence. Much less is known, however, about the influence of dynamical effects on the viscosity in its supercooled state. Here we address this issue for the first time, using molecular dynamics simulations to investigate the shear-rate dependence of the viscosity of supercooled water as described by the TIP4P/Ice model. We show the existence of a distinct cross-over from Newtonian to non-Newtonian behavior characterized by a power-law shear-thinning regime. The viscosity reduction is due to the decrease in the connectivity of the hydrogen-bond network. Moreover, the shear thinning intensifies as the degree of supercooling increases, whereas the cross-over flow rate is approximately inversely proportional to the Newtonian viscosity. These results stimulate further investigation into possible fundamental relations between these nonequilibrium effects and the quasi-static Newtonian viscosity behavior of supercooled water.

The viscosity of supercooled water has been a subject of intense study, in particular with respect to its temperature dependence. Much less is known, however, about the influence of dynamical effects on the viscosity in its supercooled state. Here we address this issue for the first time, using molecular dynamics simulations to investigate the shear-rate dependence of the viscosity of supercooled water as described by the TIP4P/Ice model. We show the existence of a distinct cross-over from Newtonian to non-Newtonian behavior characterized by a power-law shear-thinning regime. The viscosity reduction is due to the decrease in the connectivity of the hydrogen-bond network. Moreover, the shear thinning intensifies as the degree of supercooling increases, whereas the cross-over flow rate is approximately inversely proportional to the Newtonian viscosity. These results stimulate further investigation into possible fundamental relations between these nonequilibrium effects and the quasi-static Newtonian viscosity behavior of supercooled water.
Supercooled liquid water has been the subject of intense investigation for decades [1,2] and continues to attract significant attention [3][4][5]. Besides the hotly debated issue concerning the possible existence of a second critical point in the supercooled regime [6][7][8], there has been a long-standing interest in the behavior of water's viscosity below the melting temperature. Particular topics of interest include the existence of a fragileto-strong transition [9,10], the relation between viscosity and molecular diffusion [11] and the effect of pressure [12].
The viscosity η of a viscous fluid is defined as the proportionality constant between the shear stress σ and the corresponding strain rateγ according to σ = ηγ [13,14]. If, for given temperature and pressure, the relation between σ andγ is linear, i.e., η is constant, the flow behavior of the fluid is said to be Newtonian [14]. Conversely, fluids for which this linearity is violated are referred to as non-Newtonian, with colloidal suspensions, many polymer melts and granular fluids as typical examples [14,15].
Many fluids display Newtonian flow behavior for sufficiently small ratesγ. Liquid water in thermodynamic equilibrium is an example, with a viscosity that is known to be constant across several orders of magnitude oḟ γ [16]. Much less is known, however, about the dynamical effects on the viscosity of water in its supercooled state. Although its magnitude is known to rise sharply as the temperature is lowered [11,17], this increase has so far only been probed for the low-rate, Newtonian limit and the question as to whether it displays a shear-rate dependence remains open.
In this Letter we consider this issue for the first time, investigating the influence of the flow-rate on the shear viscosity of supercooled water using atomistic-level simulations. In particular, we employ non-equilibrium molec-ular dynamics (NEMD) simulations in which we impose shear deformations at a constant rateγ and measure the associated shear stress σ. To describe the interactions between the water molecules we employ the TIP4P/Ice water model [18], which is among the best molecular models for water [19] and has a melting point T m = 271 K that is close to the experimental value. All simulations have been carried out using the LAMMPS package [20]. The long-range intermolecular electrostatic interactions for the TIP4P/Ice model are calculated using the particleparticle particle-mesh (PPPM) scheme [21] and the intramolecular bond lengths and angles are held fixed using the SHAKE algorithm [22].
All the flow simulations are carried out using a computational cell containing 10800 water molecules. The cells are first allowed to equilibrate at zero external pressure and constant temperature, allowing fully flexible cells. This is achieved using a Parrinello-Rahman-type barostat [23] and a Langevin thermostat [24] with damping constants of 2 and 0.2 ps, respectively. The corresponding equations of motion are integrated using velocity-Verlet algorithm with a time step of ∆t = 1 fs. Subsequently, the nonequilibrium flow simulations are carried out at constant volume and isothermally, with temperature control implemented using a Langevin thermostat with a damping constant of 0.2 ps. The pure shear deformations are imposed using LAMMPS's fix deform command with the remap x option, allowing the molecules to adjust to the cell deformation without requiring an explicit velocity profile. This approach has shown to give good agreement with the alternative SLLOD approach [25]. Due to the appreciable cell distortions during the NEMD simulations, the reciprocal space part of the PPPM scheme is reset several times during a run, approximately after every ∼ 1% of deformation. function of the accumulated strain, γ =γ t, along six flow simulations at the deeply supercooled condition at T = 226 K.
The stress-strain curves display non-monotonic behavior that is typical of viscoelastic fluids, as has been observed in a variety of systems, both experimentally as well as in simulations [26][27][28][29][30][31][32]. At the early stages of the flow process the stress increases linearly with strain, typifying a solid-like elastic response characterized by a modulus that is independent of the deformation rate. Subsequently, the contribution of viscous relaxation processes becomes significant, first reducing the elastic increase of the shear stress to reach a maximum, σ max , followed by a final decay to a steady-state plateau value, σ ∞ . Both σ max and σ ∞ decrease as the flow rate is reduced, as the stress relaxation processes are active during longer periods of time for a given state of deformation. Indeed, foṙ γ = 2 × 10 7 s −1 the stress maximum has disappeared altogether and the stress-strain curve rises monotonically to its steady-state value.
The plateau value σ ∞ is the shear stress that is required to maintain steady-state flow at a prescribed ratė γ and the corresponding steady-state shear viscosity is then given by Fig. 2a) displays this viscosity as a function of flow rate for supercooled TIP4P/Ice water at 226 K, 246 K and 266 K, respectively. For all three temperatures the flow response can be classified into two regimes. For low rates the viscosity is independent ofγ, meaning that flow is Newtonian under these conditions. Subsequently, there is a cross-over into a non-Newtonian regime in which the viscosity decreases with growing flow rates, also known as shear thinning. Furthermore, this cross-over depends strongly on the temperature: while at 226 K non-Newtonian behavior sets in forγ 10 7 s −1 , the Newtonian flow regime persists up to flow rates ofγ ∼ 10 10 s −1 at 266 K.
To quantify the cross-over between Newtonian and non-Newtonian flow we analyze the simulation data in terms of the Carreau model [33][34][35][36], which provides a phenomenological description of shear thinning that has shown to be accurate for fluids with relatively low Newtonian viscosities, η N 1 Pa·s [36], which is the case for the present TIP4P/Ice simulations. The Carreau model treats shear flow as a stress-assisted thermally activated process involving a broad distribution of energy barriers and gives a shear viscosity that depends on the flow rate according to [33,34,36] where η N is the Newtonian viscosity,γ 0 is a characteristic cross-over rate and n is the shear-thinning exponent with a value between 0 and 1. In the limit of large flow rates this model gives rise to a power-law decay of the viscosity according to η ∞ ∼γ n−1 .
The lines in Fig. 2a) depict the least-squares regression results for the Carreau model of Eq. (2) with respect to the NEMD viscosity data. The agreement between model and simulation is very good across the entire range of flow rates for all three temperatures, clearly showing a power-law dependence of the viscosity in the shear thinning regime. The accuracy of the Carreau model can be further verified by comparing its estimate for the Newtonian viscosity η N to results from independent equilibrium calculations. Specifically, since η N represents the shear viscosity in the limit of vanishing flow rate, it can be computed using the Green-Kubo (GK) formalism [37][38][39], which expresses it in terms of stressstress autocorrelation functions that can be computed using equilibrium MD simulations. The equilibrium runs used to compute the GK viscosities are based on a cubic cell containing 2000 water molecules that are first equilibrated at zero pressure and constant temperature using the same approach used for the 10800-molecule cells. Subsequently, five independent NVT equilibrium runs are carried out to sample the components of the stress tensor and determine the stress-stress autocorrelation functions P αβ (0)P αβ (t) , where P αβ is an off-diagonal component of the stress tensor. The Green-Kubo viscosities are then computed as with V the volume of the system, T the temperature, and k B Boltzmann's constant. Aside from the three offdiagonal components P xy , P xz , and P zy , there are two other independent components, 1 2 (P xx −P yy ) and 1 2 (P yy − P zz ), that can be used due to rotational invariance [40]. Accordingly, η N is estimated using the average over these five components and over five independent equilibrium runs. Fig. 2b) presents a comparison between the NEMD Carreau results and equilibrium GK shear viscosities. The agreement is excellent for all three temperatures, providing further validation of the Carreau model as an adequate descriptor of the rate dependence of the shear viscosity in supercooled TIP4P/Ice water. A further observation based on the results in Fig. 2b) is that supercooled TIP4P/Ice water behaves as a fragile liquid for the considered temperatures [41], given that the logarithm of η N as a function of the inverse temperature 1/T is supralinear, constituting super-Arrhenius behavior.
The two other parameters of the Carreau model quantify the nature of the Newtonian to non-Newtonian transition and their behavior is plotted in Figs. 2c) and d). Fig. 2c) plots the characteristic rateγ 0 as a function of the Newtonian viscosity η N . As noted before, the transition to the shear-thinning regime sets in for lower flow rates as the temperature reduces and the Newtonian vis-cosity grows. More interestingly, the functional dependence is well described by a power law with exponent −1.16 ± 0.01, implying a direct relationship between the nonequilibrium parameterγ 0 and the equilibrium property η N . Fig. 3c) shows that the shear thinning exponent n decreases substantially as the degree of supercooling is enhanced, implying that the shear-thinning effect becomes more pronounced as the temperature is reduced. We will further discuss this point below.
There are a number of microscopic processes that can lead to the power-law viscosity behavior of the Carreau model seen in Fig. 2a) [36]. A common mechanism concerns a change in some order parameter that describes correlations between neighboring molecules [36,42]. For instance, for shear thinning in fluids composed of chain molecules, a relevant order parameter is one that measures their alignment along the flow direction [43,44]. Here, we investigate the evolution of the hydrogen bonding during the flow simulations. To determine the hydrogen-bond statistics, we adopt the definition that a HB is present whenever the distance between a proton and an oxygen satisfies 1.1Å< d OH < 2Å. Fig. 3a) displays the mean number of hydrogen bonds (HBs) per molecule, n hb , as a function of strain at T = 246 K for the flow ratesγ = 2 × 10 8 , 2.5 × 10 9 and 5 × 10 10 s −1 . These particular three values correspond to the Newtonian, the cross-over and shear-thinning regimes for this temperature, respectively. In the Newtonian regime n hb remains constant throughout the entire simulation and the connectivity of the HB network remains unaffected by the flow. As the rate increases to the Carreau crossover value, however, the steady-state HB connectivity becomes discernibly lower, reducing even further for the highest flow rate. As mentioned above, molecular alignment during the shearing process may also possibly play a role in the shear thinning, as is the case in systems where elongated molecules are involved [43,44]. To verify this possibility for water we analyze the statistics of HB directions during the shearing process. As seen in the inset of Fig. 3b), the HB direction cosines with respect to the x, y and z directions are uniformly distributed, indicating that the HB directionality is isotropic, displaying no preferred alignment direction.
These results indicate that the shear thinning arises from the reduction of HB connectivity, which is consistent with theoretical arguments [45]. The origin of this decrease and its dependence on the flow rate is associated with time-scale differences between the imposed flow and molecular rearrangements. In the Newtonian regime the latter is sufficiently short for the molecular rearrangements to accompany the imposed flow and maintain the average connectivity of the HB network. In the non-Newtonian shear-thinning regime this is no longer the case, with the molecular orientations systematically lagging behind the imposed flow, leading to the reduction of the HB connectivity in the steady state flow. This is illustrated in Fig. 3b) which depicts the time evolution of n hb along a simulation in which the system is first subjected to a constant flow rate ofγ = 5×10 10 s −1 at 246 K until reaching its steady state, after which the deformation is halted and the system is allowed to relax at a fixed cell geometry. During the flow stage the mean number of HBs per molecule rapidly decreases to its steady-state value. Subsequently, after halting the deformation, n hb relaxes to its equilibrium value by an approximately exponential process with a time constant τ m 0.07 ns. While this time scale is ∼ 300 times shorter than that associated with the lowest flow rate in Fig. 3a), it is ∼ 3 times larger compared to that of the highest.
Finally, the increasing intensity of the shear-thinning effect with reducing temperature, as reflected by the decrease of the Carreau exponent n in Fig. 2d), also correlates with the evolution of the average number of n hb . This is shown in Fig. 3c) which depicts the steady-state flow values of n hb , normalized by their equilibrium values n 0 hb , as a function of the flow rate for T = 226, 246 and 266 K, respectively. Due to the shear thinning effect, as seen in Fig. 2a), n hb decreases as the flow rate grows. Moreover, this decrease is stronger in relative terms as the temperature is lowered: whereas forγ = 5 × 10 10 s −1 a reduction of ∼ 15% with respect to its equilibrium value is observed at 226 K, it is only ∼ 6.5% at 266 K.
In conclusion, we have performed a series of NEMD simulations to investigate the shear-rate dependence of the viscosity of supercooled water as described by the TIP4P/Ice model for three different degrees of supercooling. In all cases we find a distinct Newtonian-to-shearthinning crossover that is well-described by the Carreau model. The shear-thinning effect becomes stronger as the temperature is reduced, with a thinning exponent that decreases and with non-Newtonian behavior setting in for lower deformation rates. Interestingly, the results suggest a power-law relationship between the nonequilibrium cross-over rate parameterγ 0 and the equilibrium Newtonian viscosity property η N . On the molecular scale the shear thinning correlates with a significant reduction in the connectivity of the HB network, which is associated with time-scale differences between the deformation protocol and molecular rearrangements. Moreover, the connectivity reduction increases in relative terms as the temperature is lowered, giving rise to the stronger shearthinning effect at lower temperatures.