Anomalies in supercooled NaCl aqueous solutions: A microscopic perspective

In this work we studied the effect of NaCl on the thermodynamic and dynamic properties of supercooled water, for salt concentrations between 0.19 and 1.33 mol kg(-1), using molecular dynamic simulations for TIP5P∕E water model and ion parameters specially designed to be used in combination with this potential. We studied the isobaric heat capacity (C(p)) temperature dependence and observed a maximum in C(p), occurring at T(m), that moves to lower temperature values with increasing salt concentration. Many characteristic changes were observed at scaled temperature T∕T(m) ∼ 0.96, namely a minimum in the density of the system, a reduction of the slope of the number of hydrogen bonds vs. temperature, and a crossover from Vogel-Tamman-Fulcher to Arrhenius dynamics. Finally, at low temperatures we observed that water dynamics become heterogeneous with an apparently common relationship between the fraction of immobile molecules and T/T(m) for all studied systems.


I. INTRODUCTION
In spite of being one of the most abundant substances in earth, the thermodynamic properties of water are still a matter of controversy for the scientific community.Several anomalous properties, such as a maximum density at 277 K, 1 an increase upon cooling of both the constant pressure heat capacity (C p ) (Refs.2-4) and the isothermal compressibility (κ), 5 characterize this ubiquitous substance.At normal pressure, supercooled water spontaneously freezes as it reaches the homogeneous nucleation temperature, T h = 235 K. 6 However, at very high cooling rates freezing can be avoided and supercooled water becomes a glass.Amorphous solid water can also be obtained by vapor deposition over a cold plate.The glass transition temperature, determined by slowly heating the amorphous phase, is T g ≈ 136 K. 7 The experimental determination of thermodynamic and dynamic properties of water between T h and T g is impossible in bulk samples.0][11][12][13][14] Besides, much has been known of the properties of supercooled water in this temperature region from molecular dynamic simulations [15][16][17][18] exploiting the advantage that ice nucleation is not attained in simulation time windows.
Computer simulation results show that C p (Ref.16) and κ (Ref.19) reach maximum values in the supercooled state, which move to lower temperatures with increasing pressure.a) Author to whom correspondence should be addressed.Electronic mail: longinot@qi.fcen.uba.ar.
The temperature values of maximum C p and κ define the "Widom line."Kumar et al. 16 observed that the locus of extreme C p is intimately related to the breakdown of the Stokes-Einstein relation (diffusion-viscosity correlation) and the setup of heterogeneous dynamics.Experimental results also show a crossover temperature where a change from a non-Arrhenius to an Arrhenius behavior 12,13,20 was observed.
Considering the abundance of ionic solutes in aqueous media, the properties of aqueous ionic solutions are of great interest for diverse fundamental and applied reasons.Aqueous ionic solutions at ambient temperature have been extensively studied by experimental and computational methods during many decades.However, the experimental and computational information of supercooled aqueous ionic solutions is rather scarce.
Archer and Carter 21 measured the heat capacities of pure water and NaCl aqueous solutions from 0.05 mol kg −1 to 6 mol kg −1 , at 0.1 MPa, at temperatures down to 236 K for water and 202 K for concentrated NaCl aqueous solutions.They observed, for pure water and NaCl aqueous solutions at concentrations smaller than 2 mol kg −1 , that C p increases in the supercooled region when temperature is decreased.This anomalous increment moves to lower temperature values with increasing salt concentration and disappears for salt concentrations higher than 2 mol kg −1 .It is important to note that no maximum C p is observed in the experiments because of the freezing of the system.Simulations allow probing an extended temperature range.
Recently, Corradini et al. [22][23][24] performed molecular dynamics simulations of water and NaCl solutions using the TIP4P water model.The results indicate that the presence of ions does not modify the mechanical stability limit, defined as the locus of diverging isothermal compressibility, but causes a shift of the temperature of maximum density to lower values as salt concentration is increased.
Holzmann et al. 25 have analyzed the effect of pressure and NaCl in the thermodynamic and dynamic properties of TIP4P/E water in a temperature interval that goes from 340 K down to 230 K, that is, 12 K below the melting temperature of the model.They observed that with increasing pressure and salt concentration the temperature of maximum density shifts to lower values.The results also show that application of pressure or addition of NaCl produce a reduction in the diffusion coefficient of water at high temperatures, while the opposite effect is observed at low temperatures.The authors claim that the effect of pressure and the addition of salt can be considered as "two sides of the same coin." In this work we have studied the thermodynamic and dynamic properties of water and NaCl aqueous solutions in a temperature range that includes the supercooled state.We will show that the maximum in C p observed in pure water is shifted toward a lower temperature with the addition of NaCl.Using the temperature of maximal C p to define a scaled temperature T R , given by the ratio between the temperature and the temperature of maximum C p (T R = T/T m ), we have identified several features occurring at T R ∼ 0.96 for pure water and the salt solutions.For example, a minimum in the density of the system, a reduction of the slope of the relation between the number of hydrogen bonds and temperature, and a crossover from Vogel-Tamman-Fulcher (VTF) to Arrhenius dynamics all occur for T R ∼ 0.96.Finally, at low temperature we find an heterogeneous behavior for the diffusion of water molecules that appears to be a common function of T R .
In Sec.II we describe the models and simulations details.In Sec.III we show and discuss our results.Section IV summarizes our main conclusions.

II. MODELS AND METHODS
We performed molecular dynamic simulations of neat water and NaCl aqueous solutions.All simulations were performed using GROMACS 4.04 (Ref.26) with periodic boundary conditions.The simulation box for the pure water system contained 878 water molecules.The NaCl aqueous solutions (0.19, 0.65, and 1.33 mol kg −1 ) were prepared by replacing n water molecules by n/2 Na + and n/2 Cl − ions, with n = 6, 20, and 40 for increasing salt concentration.Simulations were performed under NPT conditions with p = 0.1 MPa and temperatures between 217 and 298 K.
The water molecules were modeled using the TIP5P/E water model 27 and the ion parameters were taken from Gladich et al. 28 These parameters were designed to be used in combination with the TIP5P/E water model and carefully fitted to reproduce as best as possible experimental information for the position of the first peak of the ion-oxygen radial distribution function, the coordination number of the first solvation shell, and the solvation free energy.The details of TIP5P/E model can be found in previous publications. 27,29 he Lennard-Jones parameters for the ions are σ Na = 0.385 nm, ε Na = 2.2 J mol −1 , σ Cl = 0.385 nm, and ε Cl = 2.972 kJ mol −1 .The oxygen-ion Lennard-Jones parameters are obtained by the rules A spherical cutoff at 0.9 nm was imposed for the Lennard-Jones and short-range Coulombic interactions.The long-range electrostatic interactions were accounted for using the particle mesh Ewald method.The temperature was controlled with a Nosé-Hoover thermostat with a characteristic period of temperature fluctuations equal to 0.1 ps.Pressure was scaled using a Parinello-Rahman barostat with a reference compressibility of 4.5 × 10 −4 MPa −1 and a characteristic period for pressure fluctuations equal to 0.5 ps.Replica exchange MD simulations were performed in order to reduce the equilibration time and ensure a complete exploration of the configurational space.The typical time for a configurational swap attempt between neighboring replicas was 500 ps.All productions runs are of 40 ns or longer.
The enthalpy is determined through H = E + P V where E and V are the instantaneous total energy and volume.To calculate the diffusion coefficient, performed with trajectories coupled to a thermostat, we use r 2 = 6Dt where r 2 is the mean square displacement (msd).Swaps between neighboring replicas were taken into account in the calculation of the msd of water molecules.The msd of water molecules is calculated in 500 ps intervals discarding 10% of the data after the exchange.No corrections for system size dependence of D, as performed by Yeh and Hummer, 30 were done in this work.Diffusion constants of water were determined for all and "free" water molecules.In the latter case diffusion of molecules in the first hydration shell of the ions were discarded from the analysis.As hydration molecules we considered molecules with Na + -oxygen and Cl − -hydrogen distances smaller than the first minimum in the Na + -oxygen or Cl − -hydrogen radial distribution functions.

III. RESULTS AND DISCUSSION
We determined the temperature dependence of the constant pressure heat capacity, C p, for pure water and NaCl aqueous solutions at concentrations 0.19, 0.65, and 1.33 mol kg −1 .The calculation of C p was carried out following two different procedures: from the fluctuations of the enthalpy, and by fitting the variation of the enthalpy with temperature with a fourth order polynomial and direct differentiation of the fitting function.Both procedures gave essentially the same results; however the polynomial fit exhibits smaller scatter in C p .Therefore the values shown in Fig. 1 correspond to the second method.For all the studied systems, the C p shows a maximum in the supercooled regime.For pure water, the maximum occurs at 247 K. Experiments, limited to temperatures larger than T h , show that the C p increases with decreasing temperature without reaching a maximum.The presence of the maximum in the simulations suggests that the temperature of homogeneous nucleation for TIP5P/E must be larger than 247 K.The maximum in the heat capacity of water has also been predicted with the scaled parametric equation of  state of Fuentevilla and Anisimov, 31 who using the experimental C p values to fit its parameters, found the maximum of C p at ∼232 K. Other molecular dynamic simulations for TIP5P and ST2 potentials 16 have also seen the non-monotonic dependence of C p on temperature for pure water.
The temperatures, at which C p goes through a maximum, T m , are summarized in Table I for all NaCl aqueous solutions studied in this work.Our results show that an increment in the concentration of NaCl yields lower C p values and shifts its maximum to lower temperatures.The variation of the heat capacity curves with salt concentration are in line with the experimental observations by Archer and Carter. 21As in the case of pure water, the salt solutions crystallize before a maximum in C p could be observed, and the experiment shows a shifting of the C p vs. T curves toward lower temperatures as the salt concentration is increased.
Figure 2 shows the density of water and NaCl aqueous solutions as a function of a scaled temperature, T R , defined as T/T m .The lines correspond to fourth order polynomial fits of the raw data.Interestingly, for all the systems the density has a maximum at T R ≈ 1.13 implying that both, T m and the temperature of maximum density, have the same dependency with salt concentration.A minimum in the density at low temperature is also suggested by the simulation data for pure water.This minimum has been mentioned previously in some experimental 10,32 studies of nanoconfined water and computational studies 15,17,18 of bulk pure water.For the NaCl solutions a minimum in the density is insinuated by the simulations, except for the largest concentration that shows a continuous decrease of the density with decreasing ) FIG. 2. Density of water and NaCl aqueous solutions as a function of T R .Lines are fourth order polynomial fits of raw data, given by symbols.Concentrations in colors as in Figure 1.
temperature in the deep supercooled regime for the simulated temperature range.
Figure 3 shows the thermal expansion coefficient (α p ) of water and NaCl solutions as a function of T R , calculated from the polynomial fit of the density-temperature data.It is interesting to notice that a common behavior could be observed for all systems.The temperature of maximal density (high temperature zero of α p ) is observed at T R ≈ 1.13 for water and all the studied NaCl solutions.The low temperature zero of α p (minimum density) is observed for T R < 0.96 but it seems to depend on salt concentration.
In order to correlate the effect of salt in the thermodynamic properties with changes in the structure of water we plotted the number of hydrogen bonds per water molecule (n HB ) as a function of T R for pure water and the three studied NaCl aqueous solutions.Hydrogen bonds are defined by two commonly defined conditions.First, the oxygen-oxygen distance is within a cutoff set equal to the first minimum in the O-O radial distribution function.Second, the angle formed between the hydrogen, the donor oxygen and the acceptor oxygen is smaller than a cutoff, taken in this work as 30 • or 15 • .Figure 4 shows the dependence with temperature of the number of hydrogen bonds for a cutoff angle of 30 ) FIG. 3. Thermal expansion coefficient (α p ) of water and aqueous NaCl solutions as a function of T R .The color scheme is as in Figure 1.The corresponding figure for a cutoff angle of 15 • shows similar temperature dependence.
In Fig. 4 it can be observed that the slope of the hydrogen bonds temperature dependence decreases for T R < 0.96 for water and the two lower concentrations of NaCl studied in this work.The temperature at which the slope changes seems to be similar to that where the minimal density is predicted for the same three conditions.Interestingly the higher salt concentration does neither show the change in slope in n HB nor a minimum in the density.It should be noticed that the number of hydrogen bonds decreases with NaCl concentration since all molecules, including those in the hydration shell of the ions, were taken into account in the calculation.
It is tempting to suggest that at the scaled temperature where the variation of the number of hydrogen bonds with temperature changes slope, the solutions undergo a structural transition from a high density liquid to a more ordered low density liquid.4][35] The simulation results however, are not conclusive.Corradini et al. 36 have recently studied the effect of NaCl 0.67 mol kg −1 in the liquid-liquid transition critical point (LLCP) observed in supercooled water using molecular dynamics simulations of TIP4P water.Results show that the LLCP for NaCl aqueous solution is observed at higher temperature and lower pressure values than that observed for pure water.The authors conclude, in line with our results, that the ions stabilize the high density liquid phase shrinking in pressure the low density phase.We now turn to the changes in the dynamical properties of water in the different solutions by looking at the temperature dependence of the diffusion constant (D). Figure 5 displays the calculated values of D in Arrhenius form, for absolute (a) and scaled temperatures (b).The diffusion coefficient temperature dependence is non-Arrhenius and, in all but the lowest temperatures, can be described by the VTF equation, [37][38][39] where D 0 and C are fitting parameters and T 0 corresponds to the temperature at which dynamics would be completely arrested.However, below T 0 dynamics deviate from the VTF behavior.Fits of the diffusion constants to the open symbols in Figure 5(a) with Eq. ( 2) are displayed as full lines in the figure and the fitting parameters are shown in Table II.Results show that T 0 decreases and C increases with increasing NaCl, that is, complete arrest of the dynamics, according to VTF, takes place at lower temperatures with increasing NaCl concentration.
For water, at 1000/T > 4.2 K −1 (T < 238 K) the diffusion deviates from VTF behavior.It is tempting to describe the behavior of D at T below 238 K as Arrhenius.The fit is shown in Figure 5(a) as the dashed line with an activation energy of 72 ± 3 kJ mol −1 .Mallamace et al. 12 measured the diffusion coefficient of water in confined environments and observed a crossover from VTF to Arrhenius and in the latter found an activation energy of 20 kJ mol −1 .TIP5P/E gives diffusion coefficients that decrease faster with temperature than measured data, 40,41 and that behavior may explain the activation energy difference with the experimental data.Experimental data for supercooled water in confined systems reported by Liu et al. 9 show a change from a VTF to an Arrhenius behavior for the translational relaxation time at 224 K at ambient pressure.Simulation results for SPC/E water in confined silica 42 show the same transition at 215 K, which coincides with a peak in the specific heat capacity at the same temperature, in line with our results.
In Fig. 5(b) we plotted the diffusion as a function of the inverse scaled temperature, 1/T R , and it can be seen that in the non-VTF regime all the curves collapse into a universal function.This result suggests that the water dynamics are Arrhenius-like, with a universal (scaled and dimensionless) activation energy given by 72 000/(R T m 0 ), where T m 0 is the temperature of maximum C p of pure water and R is the gas constant given in J mol −1 .
Figure 6 shows the ratio between the diffusion coefficients of water in NaCl solutions and in pure water, D/D w , as a function of salt concentration at several temperatures.Above the freezing temperature, the diffusion constant decreases with increasing NaCl concentration, while below this temperature the opposite behavior is observed.These predictions were also observed for TIP5P water by Kim and Yethiraj 43 over a more limited range of temperatures and an extended concentration range.Holtzmann et al. 25 have also reported the same crossover for TIP4P/E potential around the freezing temperature, but shifted to lower temperature since the freezing point for TIP4P/E is lower than that for  TIP5P model.This behavior has been explained by a change of the effect of the salt from a structure-making character at room temperature, to a structure-breaking character at low temperatures. 25,43  the simulations the changes in the diffusion of water with temperature seem to arise from two effects.One is the natural decrease of mobility with temperature and the other is the fraction of water molecules that become immobile throughout the length of the simulation.Figure 7 shows the fraction of immobile water molecules (f 0 ) as a function of the scaled temperature.As immobile molecules we considered those water molecules whose diffusion constants were smaller than 10 −11 cm 2 s −1 during the whole simulation time (at least 40 ns for all simulations).Above T R = 1 all water molecules are mobile.For T R < 1 the fraction of immobile molecules increases with decreasing temperature, and within the uncertainties of the simulations it shows common trend.This common dependence suggests that the water molecules in the first hydration shell of the ions behave in a similar dynamical way than the rest of the molecules.To test this hypothesis we calculated the diffusion coefficient for the water molecules that do not participate on the first hydration shell during the whole simulation time.Results show that the temperature dependence of the diffusion constant of "free" molecules is similar to that observed for all molecules.The crossover temperature where the diffusion deviates from the VTF behavior and where the dynamics become heterogeneous are the same in both cases.It should be stressed that for NaCl 0.65 mol kg −1 between 235 and 225 K around 18% of all molecules participate in the first hydration shell of the ions, while for NaCl 1.33 mol kg −1 this number increases to 34% and 46%, respectively for the same mentioned temperatures.Therefore, although the fraction of molecules in the hydration shell of the ions is quite big, the difference between the temperature dependence of the diffusion constant of all molecules and "free" molecules is negligible.This confirms our hypothesis that molecules in the first hydration shell of the ions should behave in a similar dynamical way than the rest of the molecules.

IV. CONCLUSIONS
We performed an exhaustive analysis of the thermodynamic, structure, and dynamic properties of supercooled water and NaCl aqueous solutions of varying concentration (0.19, 0.65, and 1.33 mol kg −1 ).Our work is based in molecular dynamic simulations of TIP5P/E water in the NPT ensemble at atmospheric pressure.
We observed the maximum in C p for water in the supercooled regime, in line with previous simulation works.We studied the effect of NaCl on the heat capacity and showed that the maximum in C p shifts to lower temperatures with increasing salt concentration.In order to analyze the results, we found it convenient to introduce the scaled temperature T R , defined as the temperature scaled by the temperature of maximum C p .The temperature of maximum density occurs at T R ∼ 1.13 for pure water and all the studied solutions.A minimum in the density is observed for T R < 0.96 and the slope of the temperature dependence of the number of hydrogen bonds decreases for T R < 0.96 for water and the two more dilute solutions (NaCl concentrations 0.19 and 0.65 mol kg −1 ).These findings suggest that a structural broad transition between a high density liquid and a low density liquid, already proposed by several authors, occurs around 0.96 T R .
The diffusion coefficient of water, both in the pure liquid and in the salt solutions, shows a crossover from VTF to Arrhenius behavior as the temperature is decreased.Interestingly, the Arrhenius behavior appears as a universal function in terms of the scaled temperature T R .Moreover, the scaled transition temperature is T R ∼ 0.96, a temperature that characterizes some of the static features discussed above.The mobility of the molecules that do not participate in the ions first hydration shell is essentially the same as that of the total system, for all studied temperatures.We also observed a heterogeneous dynamics behavior of the water molecules at low temperatures.For T R < 1, a fraction of water molecules remain immobile during the whole simulation.This fraction of immobile molecules follows an apparently common behavior, independent of salt concentration.This result suggests that the water molecules in the first hydration shell of the ions have a similar dynamics behavior than the rest of the system.

FIG. 4 .
FIG. 4. Number of hydrogen bonds per water molecule (n HB ) for H 2 O andNaCl aqueous solutions as a function of T R for a cutoff angle of 30 • .The color scheme is the same as in Figure1.

FIG. 5 .
FIG.5.Arrhenius plot of the diffusion coefficient of water in pure water and NaCl aqueous solutions as a function of temperature (a) and scaled temperature, T R (b).Colors are the same as in Figure1.(a) Symbols correspond to raw data and full lines to fitting curves given by Eq. (2) for data plotted with empty symbols.Full symbols were discarded from the fit since deviations to Eq. (2) are observed for those data.For water the dashed line represents the fit of full symbols with the Arrhenius equation.(b) Lines are given as a guide to the eye.

FIG. 7 .
FIG.7.Fraction of immobile water molecules (f 0 ) as a function of the scaled temperature for water and NaCl aqueous solutions.Symbols correspond to raw data and lines are used as a guide to the eye.Colors are the same as in Figure1.

TABLE I .
Temperature of maximum C p (T m ) and maximum C p values (C p,m ) for water and NaCl aqueous solutions.

TABLE II .
Parameters of the VTF equation (Eq.(2)) for pure water and NaCl aqueous solutions.