Salt aqueous solutions are relevant in many fields, ranging from biological systems to seawater. Thus, the availability of a force-field that is able to reproduce the thermodynamic and dynamic behavior of salt aqueous solutions would be of great interest. Unfortunately, this has been proven challenging, and most of the existing force-fields fail to reproduce much of their behavior. In particular, the diffusion of water or the salt solubility are often not well reproduced by most of the existing force-fields. Recently, the Madrid-2019 model was proposed, and it was shown that this force-field, which uses the TIP4P/2005 model for water and non-integer charges for the ions, provides a good description of a large number of properties, including the solution densities, viscosities, and the diffusion of water. In this work, we assess the performance of this force-field on the evaluation of the freezing point depression. Although the freezing point depression is a colligative property that at low salt concentrations depends solely on properties of pure water, a good model for the electrolytes is needed to accurately predict the freezing point depression at moderate and high salt concentrations. The coexistence line between ice and several salt aqueous solutions (NaCl, KCl, LiCl, MgCl2, and Li2SO4) up to the eutectic point is estimated from direct coexistence molecular dynamics simulations. Our results show that this force-field reproduces fairly well the experimentally measured freezing point depression with respect to pure water freezing for all the salts and at all the compositions considered.

When a solute is added to water, the freezing point of water decreases. The difference between the freezing temperature of pure water and that of a solution is denoted as freezing point depression. The freezing point depression caused by an electrolyte is larger than that of a non-electrolyte at the same molality (i.e., moles of solute per kilogram of solvent). This is due to the fact that electrolytes when dissolved into water dissociate into ions, thus generating more particles in solution. If ideal behavior is considered, the freezing point depression is proportional to the number of ions generated when dissolved (i.e., 2 for a 1:1 electrolyte and 3 for a 1:2 electrolyte). The constant of proportionality is related to the properties of pure water and its experimental value is 1.86 K per molal unit. Not surprisingly, electrolytes are, in fact, often used to avoid the formation of ice in roads, since when salt is added, water only freezes at temperatures below 0 °C. In the particular case of NaCl, the largest freezing point depression is reached at the eutectic point, which is found at about 21 K below the melting point of pure water. This is the minimum temperature at which a NaCl aqueous solution could be present as the stable phase. At room pressure, the solubility of salts in ice is practically zero, and for this reason, when water freezes from an ionic solution, one obtains almost pure ice (a quite small amount of ions can incorporate into the ice at the micromolar range at best1,2). At the freezing temperature, the chemical potential of water in pure ice is equal to that of water in the ionic solution. Thus, the freezing point depression gives some information about the activity of water in the ionic solution.

Liquids can be supercooled,3 and this is true for pure water and for ionic solutions. When liquids are broken in small droplets of the order of a few microns of radius, they can be supercooled up to 40 K below the freezing temperature. The degree of supercooling depends on the size of the droplets and on the velocity of cooling. Sometimes, the temperature at which these droplets freeze is denoted as the freezing temperature (although this is a kinetic definition, not a thermodynamic one, and for this reason, the term homogeneous nucleation temperature is preferred). Nevertheless, the term melting temperature is reserved for the temperature at which the chemical potential of water is identical in the two phases (ice and solution).4 In this work, we shall use the concept freezing temperature in a thermodynamic way (i.e., where the chemical potential becomes identical in both phases), which corresponds to that denoted as melting temperature in kinetic studies of freezing.

Determining the freezing point depression from computer simulations is challenging from a technical point of view. A possible route is to determine the chemical potential of ice (for a certain pressure) as a function of the temperature and that of water in the solution (for a certain pressure and composition) as a function of the temperature and to analyze the temperature at which both chemical potentials are identical. This route is certainly possible but somewhat cumbersome. One needs to evaluate the chemical potential of water in the solid phase using, for instance, Einstein crystal methods5,6 and the chemical potential of water in solution. Determining the chemical potential of water in an ionic solution (and the chemical potential of the salt) is also possible as it has been shown over the last few years, but it is certainly a challenging and rather difficult calculation.7,8 One wonders if another simpler route would be possible. An interesting possibility is the use of the direct coexistence method.9,10 This method has proved to be successful in determining the freezing temperature of pure water, and one could imagine that the method could be applied to electrolyte solutions. In fact, Jungwirth started fifteen years ago to study the freezing of water from an ionic solution1,11 (paying special attention to brine rejection from the ice), and this has been followed by other authors.12–14 In 2008, Kim and Yethiraj15 performed a pioneering study where the freezing point depression was determined from computer simulations using the direct coexistence method, putting together ice and a concentrated solution of the electrolyte and waiting to reach the equilibrium running for about 50 ns (a tour of force at that time). More recently, Conde, Rovere, and Gallo also illustrated the feasibility of direct coexistence studies to obtain the freezing point depression16 by using runs of up to 300 ns. In this case, instead of waiting to reach equilibrium, it was determined if the ice grew or melted for a certain initial concentration of salt. A similar approach was used by Soria et al.17 who also showed that the freezing temperature found from free energy calculations was identical (within statistical uncertainty) to that obtained from direct coexistence simulations for a certain concentration of salt. In the work of Soria et al.,17 kinetic aspects of the nucleation of ice in salty water were also studied, and nucleation rates were estimated using the seeding technique.18 

Although these pioneering works have illustrated the possibility of using direct coexistence as a systematic study of the relation between the concentration of the electrolyte and the freezing temperature, it has not yet been applied to different salts and types of electrolytes. Most of the force fields proposed so far for electrolytes were designed to reproduce the density of the solution and the hydration free energy using rigid non-polarizable models to describe water, for instance, TIP3P,19 TIP4P-Ew,20 or SPC/E.21 For the ions, a Lennard-Jones (LJ) center is typically used with an integer charge (in electron units, e). Probably the force field proposed by Joung and Cheatham is now the most popular.22 Properties such as solubility, viscosity, and diffusion coefficient of water were not used as target properties. How these force fields perform for these properties? The short answer is, not too well.

Concerning the solubility, it has become clear that most of the force fields predict a rather low solubility (sometimes being twice smaller than the experimental one and sometimes being ten times smaller).23–27 It has also fallen into place that the ion-pairing and precipitation often observed in simulations of ionic solutions were due to a too low solubility of the force field.28–31 In fact, the solubility of a salt in water, for a certain force field, is not necessarily equal to the experimental one.8 In general, the trend is that the solubility is too low, and the ion pairing is too high. Concerning the viscosity, it is getting across that common force fields overestimate the increase in the viscosity due to the addition of salt when compared to experiments.32,33 Finally, concerning the diffusion coefficient of water, it was pointed out by Kim et al.34 that most of the force fields (including polarizable ones) produce a larger decrease on the diffusion coefficient of water than the experimental one. How to remedy this situation? In 2014, Kann and Skinner35 suggested that using scaled charges for the ions (i.e., for a 1:1 electrolyte assigning charges smaller than one in electron units) improved the description of the diffusion coefficient of water in ionic solutions. The use of scaled charges was first advocated by Leontyev and Stuchebrukhov36–38 and was denoted as the electronic continuum correction (ECC). They suggested to use a scaled charge of 0.75 e that is obtained from a theoretical argument, which requires the dielectric constant of water at very high frequencies (when only electrons can reply to a fast oscillating electric field). This idea was also advocated by Jungwirth and co-workers39–43 and also implemented by several groups.44–47 Following this line of research in 2017, we proposed a force field for NaCl in water48 (with water described by the TIP4P/2005 model49) that used the concept of scaled charges. Our motivation was that may be the charges used to describe the dipole moment surface were not convenient to describe the potential energy surface.50,51 In 2019, we extended this idea to many different electrolytes (with stoichiometry 1:1 and 1:2), and the force field was denoted by Madrid-2019.33,52,53 The value of the scaled charge in this case was 0.85 e because this value allows us to recover the correct Debye–Huckel law, as shown by Kann and Skinner,35 i.e., it compensates the deficiencies in the dielectric constant of TIP4P/2005. It was found that the use of scaled charges improved the description of the viscosities, diffusion coefficients of water, activity coefficients, and solubilities. For example, the solubility of NaCl for the Madrid-2017 model48 is estimated to be 5.7 mol kg−1 at room conditions, in reasonable agreement with the experimental solubility, 6.1 mol kg−1.54 We expect a similar value for the Madrid-2019 model, as the change in the parameters between both models is small. For the remaining salts, the solubility has not yet been calculated. We should note that Yagasaki et al.55 recently showed that, when properly designed, force fields with integer charges can also reproduce the experimental value of the solubility. As a counterpart, the use of scaled charges deteriorates the prediction of hydration free energies56 and is not recommended either for molten salts or for the solid phase.8,57

The idea of scaled charges finds a molecular basis on the fact that the calculated charge of the ions from ab initio studies is typically in the range of 0.8 e (for anions) and 0.9 e (for cations).58,59 Thus, it seems that the use of non-integer charges (an approximation often used by the community that performs simulations of ionic liquids60) seems to improve the description of electrolyte solutions. The scaled charge approach (or ECC) can be regarded as an approximation of zero order. As a more elaborated approach, one could use models that explicitly use the concept of charge transfer between the ions and water, as it has been shown by Soniat et al.61,62 and by Yao et al.59 

In this work, we have two goals. The first one is methodological. We will explore the possibility of determining the freezing point depression by using the direct coexistence method to determine if this technique can be used to obtain reliable results using the current computational resources, which allow runs of the order of the microsecond. We shall follow the methodology of Kim and Yethiraj15 so that ice will be put in contact with an ionic solution and simulations will be performed until the system reaches the equilibrium. For a two component system, this is a quite useful approach, since according to the Gibbs phase rule, two phases can be at equilibrium for an imposed value of the pressure and temperature. That was also illustrated by Espinosa et al.27 to determine the solubility of a salt and by Kolafa63 for a system of a solid (NaCl) and an ionic solution (NaCl in water) where simulations were performed until the system reached the equilibrium.

It is interesting to point out that the melting temperature of ice Ih for the TIP3P model is around 150 K,64 and for SPC/E, it is 215 K.65 Therefore, the simulations to study the freezing point depression should be performed at extremely low temperatures, where the diffusion of water will be terribly slow and the direct coexistence method would be in trouble. Thus, it does not seem a good idea to use these water models (neither their force fields) for this purpose. A better choice would be the TIP4P-Ew model20 with a melting point of about 244 K66 and an even better choice would be the TIP4P/2005 with a melting temperature of 250 K.67 This model is convenient as water is still able to diffuse reasonably well in the range 210–250 K (i.e., for supercoolings up to 40 K). Thus, the use of TIP4P/2005 for computer simulation studies of freezing point depression is a good choice for two reasons: its melting temperature is not too low, and a force field for electrolytes (Madrid-2019) is available for this water model. TIP5P68 and/or TIP4P/Ice69 could also be considered as promising models for freezing point depression as their melting points are quite close to the experimental one. However, there are not force fields of salts explicitly developed for these water models, and in the case of TIP5P, ice II seems to be more stable64 than ice Ih at room pressure.

There is a second goal in this work. We want to test if scaled charges are able or not to predict reasonably well the freezing point depression of different salts. It seems interesting to analyze if the improvement gained by using scaled charges in properties such as viscosity, solubility and diffusion coefficient of water has as a consequence a poor performance of other properties, such as the freezing point depression. As it will be shown in this work, this fear is not justified and the model performs quite well. The good agreement found with experiments also provides an indirect indication that the activities of water should not be too far away from the experimental ones. We hope that freezing point depression will be incorporated as another test property to be checked when proposing a new force field for electrolytes in water.

The freezing point depressions of the salt aqueous solutions are estimated in this work using the direct coexistence technique.9,70 Briefly, this method consists in building a simulation box in which two phases (in this case, ice and the electrolyte aqueous solution) are initially placed next to each other and following the evolution of the system at different thermodynamic conditions. Depending on the thermodynamic ensemble used in the simulations and the number of components in the system, the fate of the interface changes: it can either reach a stable state or irreversibly evolve toward one of the two phases.

At a given thermodynamic state, two (or more) phases (I, II) are in equilibrium when the chemical potential of all the components of the system (two in our case) are equal in the coexisting phases,

μ1I(p,T,x1I)=μ1II(p,T,x1II),μ2I(p,T,x2I)=μ2II(p,T,x2II).
(1)

Here, μiα and xiα are, respectively, the chemical potential and mole fraction of component i in phase α.

The Gibbs phase rule imposes constraints on the number of phases that can coexist and applies both to single- and multiple-component systems. It states that for a system with C components and P phases, the number of independent variables or degrees of freedom F is given by

F+P=C+2.
(2)

Consequently, in one-component systems (C = 1), if both p and T are imposed (F = 2), only one phase can be stable. This explains why in an one-component system an interface cannot be stabilized in NpT simulations. In the NVT ensemble, though, only the temperature is fixed (F = 1), and thus, in these conditions, two phases of a one-component system can coexist.71 On the contrary, in a two-component system (C = 2), as, for example, the salt aqueous solutions considered in this work, we can have two phases in equilibrium even when both p and T are fixed (F = 2). In this case, the compositions of the two phases in coexistence evolve until the chemical potentials of each component are equal in the two phases [Eq. (1)]. Even though we have written the condition of equal chemical potential of component 2 in terms of its mole fraction (x2I and x2II), it could have been equally written in terms of the mole fractions of component 1 (x1I and x1II) because the mole fractions add up to one in each phase (x1I + x2I = 1 and x1II + x2II = 1). Thus, in two-component systems, there is only one independent variable that adjusts the compositions of the two coexistence phases. Note that it would also be possible to perform the direct coexistence simulations in other ensembles, such as NVE. In this case, the pressure, temperature, and composition will evolve to equilibrium. However, using the NVE ensemble has the disadvantage that special care must be taken so that the dimensions of the unit cell in the directions parallel to the interface are chosen commensurate with the equilibrium unit cell at the equilibrium temperature and pressure.71 

The fact that in two-component systems two phases can be in equilibrium at fixed p and T can be exploited in the direct coexistence method to calculate the equilibrium compositions at coexistence.27 One simply needs to build a simulation box in which the two phases are initially in contact with each other and let the system evolve. The ice slab will grow or melt depending on the initial conditions, releasing/absorbing water molecules to/from the salt solution, until the liquid phase reaches the equilibrium composition. The only subtlety is that in finite systems as those used in simulations, the initial conditions must be close enough to coexistence so that the system is able to reach equilibrium before exhaustion of any of the coexistence phases. As mentioned in the Introduction, an alternative is to follow the evolution of the interface for an initial salt concentration to bracket the concentration at equilibrium.16 In this work, we will follow the equilibrium route.

The interactions in the salt aqueous solutions were described using the recently proposed Madrid-2019 model.33 This force field describes water and the sulfate ions as rigid and non-polarizable. Water is modeled with the TIP4P/2005 potential49 and ions are represented by Lennard-Jones centers and scaled point charges (0.85 e at the monovalent and 1.7 e at the divalent ions) that, in an effective way, account for the polarization effects in the solution.

The initial configuration was built by placing a slab of 2048 molecules of ice Ih in contact with a salt aqueous solution containing 3330 water molecules and a varying number of ions depending on the temperature. The initial configuration of ice Ih with proton disorder and almost zero dipole moment satisfying the Bernal–Fowler72 rules was generated using the algorithm proposed by Buch et al.73 Initial solution concentrations between 1 and 6 mol kg−1 were considered to cover temperatures within the range T = 200–247 K. The corresponding numbers of ions (Li+, Na+, K+, Mg2+, Cl, and SO42) are given in Table I. The block of ice was placed exposing the basal plane to the salt solution and removing the overlaps between particles that lead to high repulsive forces. The freezing point does not depend on which crystallographic plane is exposed; it only changes the speed of growth/melting of the surface.74,75 The edges of the simulation box are approximately Lx ≈ 3.6 nm, Ly ≈ 3.1 nm, and Lz ≈ 15.5 nm. The ice–solution interface was aligned perpendicularly to the z axis. The dimensions of the ice slab are 3.6 × 3.1 × 5.9 nm3. Note that the width of the ice slab is large enough to avoid finite size effects on the direct coexistence method.27 

TABLE I.

Number of ions in the systems depending on the salt compound stoichiometry and the initial molality, mi. Nan is the number of anions and Ncat is the number of cations. Compounds with stoichiometry 1:1 are NaCl, KCl, and LiCl, and compound with stoichiometry 1:2 is MgCl2 and with stoichiometry 2:1 is Li2SO4.

mi (mol kg−1)Nan 1:1 and 2:1Ncat 1:1 and 1:2Nan 1:2Ncat 2:1
60 ⋯ 
120 240 
180 360 
240 ⋯ 
360 ⋯ 
mi (mol kg−1)Nan 1:1 and 2:1Ncat 1:1 and 1:2Nan 1:2Ncat 2:1
60 ⋯ 
120 240 
180 360 
240 ⋯ 
360 ⋯ 

The system is then allowed to evolve. The simulations are performed using the molecular dynamics GROMACS package76 in the NpT ensemble so that the ice slab is able to adjust its lattice constants by varying the edges of the simulation box in the two dimensions parallel to the interface. Previous works have shown that even though the pressure in the longitudinal and transversal directions to the interface is not equal due to the presence of the interface, direct coexistence simulations in the NpT ensemble give the correct results if the longitudinal dimension of the simulation cell is sufficiently large so that the contribution of the interface to the pressure tensor is small.71,77 The pressure was set to p = 1 bar and was controlled using an anisotropic Parrinello–Rahman barostat,78 with a relaxation time of 2 ps and a compressibility of 10−6. The temperature was set using the Nosé–Hoover thermostat79 with a relaxation time of 2 ps. The equations of motion were integrated using the velocity-Verlet algorithm using a time step of 2 fs. To deal with the long range electrostatic interactions, the particle mesh Ewald summations80 were used. The cut-off distance for the dispersive and the real part of the electrostatic interactions was set to 10 Å. Standard long range corrections to the energy and pressure (LRCs) were added to the calculation of Lennard-Jones (LJ) interactions. LRCs were calculated using the mean density of the system. Although this is not strictly correct for inhomogeneous systems, the densities of the solution and of ice are rather similar. The geometry of the water molecules was constrained with LINCS using a sixth order expansion and a maximum of four iterations,81,82 except for the sulfate solution for which SHAKE83 was used.

The initial interface will evolve to reach the equilibrium concentration at the simulated temperature and pressure. If the concentration of the solution is higher than that at coexistence, the ice slab will melt, releasing water to the solution to reduce its salt concentration. On the contrary, if the initial concentration is lower than that at equilibrium, the ice slab will grow to increase the concentration of the solution. The initial and final configurations of a sample interface in which an ice slab is in contact with a LiCl solution at T = 220 K are shown in Fig. 1.

FIG. 1.

Initial (top) and final (bottom) configuration of an ice slab in contact with a LiCl solution at T = 220 K and room pressure. The initial salt concentration was set to m = 4 mol kg−1, which is below the equilibrium concentration at those thermodynamic conditions. In this case, the ice slab grows, incorporating some solvent from the solution until the salt concentration of the solution reaches the coexistence value at the chosen conditions of temperature and pressure, which is roughly m = 4.48 mol kg−1.

FIG. 1.

Initial (top) and final (bottom) configuration of an ice slab in contact with a LiCl solution at T = 220 K and room pressure. The initial salt concentration was set to m = 4 mol kg−1, which is below the equilibrium concentration at those thermodynamic conditions. In this case, the ice slab grows, incorporating some solvent from the solution until the salt concentration of the solution reaches the coexistence value at the chosen conditions of temperature and pressure, which is roughly m = 4.48 mol kg−1.

Close modal

We consider that the system has reached equilibrium when the potential energy of the system and the solution salt concentration remain, on average, constant. The concentration of the solution at equilibrium was estimated from the density profiles of water [ρwat(z)] and salt [ρsalt(z)] measured along the z-coordinate that, in the chosen coordinate system, runs in the direction perpendicular to the ice–solution interface. The density profiles were evaluated by dividing the simulation box in 100 slabs of width Δz = Lz/100 and measuring the water or salt density in each of them,

ρi(z)=wi(z+Δz)LxLyΔz,
(3)

where wi(z + Δz) is the mass of the water molecules (wwat) or of the ions (wions) whose positions along the z-coordinate are between z − Δz/2 and z + Δz/2 and ⟨Lx⟩ and ⟨Ly⟩ are the average edge lengths of the simulation box in the directions parallel to the interface. These density profiles were averaged over periods of 50 ns (500 frames). As an example, the equilibrium density profile of the MgCl2 system at T = 215 K is shown in Fig. 2(a). From these plots, it is easy to determine the region in which the solution exhibits bulk behavior. Note that densities deviate significantly from the bulk at distances close to the interface, and these regions should not be taken into account in the estimates of the salt concentration in the liquid phase. From these densities, one then can easily compute the salt concentration, e.g., the molality m(z), using

m(z)=ρsalt(z)ρwat(z)×Msalt,
(4)

where Msalt is the molar mass of the ionic compound. As density profiles are measured in kg/m3, the molar mass Msalt must be expressed in kg mol−1 units. The molality of the solution, m, is then estimated by averaging m(z) over the slabs in which the solution behaves as bulk. As mentioned before, at equilibrium, the molality remains stationary with time, fluctuating about the average equilibrium value. An example of the evolution of molality with time is shown in Fig. 3. As can be seen, the salt concentration undergoes large fluctuations at short times, but oscillations become smaller and the molality adopts an average constant value at long times. The red vertical line marks the point from which the molality was averaged in the MgCl2 system at T = 215 K. This vertical line represents the time required to reach equilibrium (i.e., where fluctuations in the concentration of the solution are not larger than around 0.1 mol kg−1 of the equilibrium value). The uncertainty in the equilibrium salt concentration was estimated from the dispersion of the molality measured in these blocks. The error in the freezing point depression was estimated by assuming an uncertainty of 0.5 K in the freezing point of pure water, and adding to this an estimation of the uncertainty of the freezing point of the salt from the uncertainty in the molality. The total error in the given freezing point depression is about 0.75–1 K.

FIG. 2.

Density profiles of water and salt in the MgCl2 system at T = 215 K and an initial concentration of salt 3 mol kg−1. (a) The region comprised between the two vertical green lines was considered in the calculation of the molality. Results are averaged at times between 300 and 350 ns and the system is divided lengthwise in 100 bins. (b) Concentration profiles (in mol dm−3) of each ion are shown separately to show that the system fulfills the requirement of local electrostatic neutrality. Results are averaged at times between 1000 and 1050 ns and 30 bins are considered. For aiding visualization, the concentration of Cl is divided by 2.

FIG. 2.

Density profiles of water and salt in the MgCl2 system at T = 215 K and an initial concentration of salt 3 mol kg−1. (a) The region comprised between the two vertical green lines was considered in the calculation of the molality. Results are averaged at times between 300 and 350 ns and the system is divided lengthwise in 100 bins. (b) Concentration profiles (in mol dm−3) of each ion are shown separately to show that the system fulfills the requirement of local electrostatic neutrality. Results are averaged at times between 1000 and 1050 ns and 30 bins are considered. For aiding visualization, the concentration of Cl is divided by 2.

Close modal
FIG. 3.

Molality as a function of time, t, in the solution of MgCl2 with an initial concentration of 2 mol kg−1 in contact with an ice slab at T = 232.5 K. The red line indicates where the equilibrium is reached (i.e., when fluctuations in the concentration of the solution are not larger than around 0.1 mol kg−1 of the equilibrium value, which is represented with a blue line).

FIG. 3.

Molality as a function of time, t, in the solution of MgCl2 with an initial concentration of 2 mol kg−1 in contact with an ice slab at T = 232.5 K. The red line indicates where the equilibrium is reached (i.e., when fluctuations in the concentration of the solution are not larger than around 0.1 mol kg−1 of the equilibrium value, which is represented with a blue line).

Close modal

The profiles of the molarities of the ions (mol dm−3) are also useful are also useful to check that the system fulfills the requirement of local electrostatic neutrality, i.e., the average net charge is zero in any slab sufficiently away from the interface. Molarity profiles are shown in Fig. 2(b). Although in the proximity of the interface with ice electroneutrality is not strictly satisfied (due to the different adsorptions of cations and anions), it can be seen that electroneutrality is satisfied in the bulk solution far away from the interface.

The time needed to reach equilibrium is stochastic, and independent simulations at the same thermodynamic conditions and initial salt concentration can reach equilibrium at different times.27 However, on average, simulations starting from conditions close to coexistence will generally equilibrate faster than those starting from conditions far away from the freezing point. If the chosen initial concentrations are too far from equilibrium, the whole ice slab can melt before reaching the equilibrium composition of the salt aqueous solution, or it can grow leaving a very small region of liquid in which the bulk region is too narrow to properly evaluate the equilibrium salt concentration. Initial and final configurations of simulations exemplifying these two scenarios are shown in Fig. 4. In either of these two situations, a new simulation box must be built but tweaking the initial composition of the fluid: if all the ice slab has melted, the salt concentration of the solution must be reduced, and on the contrary, if the solution region is too narrow, the initial salt concentration of the solution has to be increased.

FIG. 4.

(a) Initial (top) and final (bottom) configurations of an ice slab in contact with a NaCl solution at T = 243 K and room pressure. The initial salt concentration was set to m = 6 mol kg−1, which is much higher than the equilibrium concentration at those thermodynamic conditions (m = 2.13 mol kg−1). In this case, the ice slab melts completely. (b) Initial (top) and final (bottom) configurations of the same system as before but now setting the initial salt concentration to m = 0.5 mol kg−1, which is much below the equilibrium concentration at those thermodynamic conditions. Now, most of the system crystallizes, leaving a fairly narrow region for the solution phase.

FIG. 4.

(a) Initial (top) and final (bottom) configurations of an ice slab in contact with a NaCl solution at T = 243 K and room pressure. The initial salt concentration was set to m = 6 mol kg−1, which is much higher than the equilibrium concentration at those thermodynamic conditions (m = 2.13 mol kg−1). In this case, the ice slab melts completely. (b) Initial (top) and final (bottom) configurations of the same system as before but now setting the initial salt concentration to m = 0.5 mol kg−1, which is much below the equilibrium concentration at those thermodynamic conditions. Now, most of the system crystallizes, leaving a fairly narrow region for the solution phase.

Close modal

When comparing with experiments, the melting temperature will be expressed as ΔT=TmsolutionTm, where Tm is the melting point of pure water. In this way, we account for the decrease in the melting point of water by the addition of salt. Note that the TIP4P/2005 exhibits a melting point of 250 K,67 i.e., about 23 K lower than the experimental value, 273.15 K. We will see that even though the absolute temperature at which freezing occurs for a given salt concentration is underestimated in absolute value, an excellent agreement is recovered when the temperature is measured with respect to the melting temperature of pure water.

Once the melting point for a given salt at a given concentration is known, as a by-product of the direct coexistence simulations, we can also obtain the water activity or, equivalently, the water activity coefficient at coexistence. The activity coefficient measures the deviation of the chemical potential of water in the solution with respect to an ideal solution.

As stated in Eq. (1), at coexistence, the chemical potentials of the two components are equal in the two phases. The chemical potential of water in the solution can be expressed in terms of the water activity (awat),

μwatsolution(T)=μwat0(T)+RTlnawat=μwat0(T)+RTln(γwatxwat),
(5)

where μwat0(T) is the chemical potential of pure water at the same conditions of p and T as the mixture. The chemical potentials are a function of the temperature and the pressure, but we leave out the pressure because all the simulations were done at 1 bar. In the last equality of Eq. (5), the activity is expressed in terms of the water molar fraction xwat and the activity coefficient γwat. At coexistence, the chemical potential of water in the solid phase is equal to the chemical potential of water in the solution [Eq. (1)]. Assuming that the chemical potential of water in the solid phase is identical to that of pure ice (given the very small solubility of salts in ice) and taking into account the condition of chemical equilibrium (i.e., the same chemical potential for water in ice and in the solution at coexistence), one obtains

μice0(T)=μwat0(T)+RTlnawat=μwat0(T)+RTln(γwatxwat).
(6)

If the mixture behaves as an ideal mixture, the activity coefficient is equal to one, and the previous expression reduces to

Δμ0(T)=μice0(T)μwat0(T)=RTlnxwat,
(7)

where the water molar fraction must be calculated taking into account the number of ions (instead of the number of salt molecules),84 

xwat=nwatnwat+nions,++nions,,
(8)

where nwat is the number of moles of water and nions,+ and nions,− are, respectively, the number of moles of positively and negatively charged ions. In Eq. (7), Δμ0(T)=μice0(T)μwat0(T) can be obtained by thermodynamic integration from the freezing point of pure water using the expression

Δμ0(T)RT=TmTHice0(T)Hwat0(T)RT2dT.
(9)

Here, Hice0 and Hwat0 are the enthalpy (per mol) of ice and liquid water in the pure system, respectively. To obtain this expression, we have used that at the melting temperature of the pure system, Tm, the chemical potentials of ice and liquid water are equal. As a first approximation, the freezing point depression of an ideal electrolyte can be computed, assuming that the integrand of Eq. (9) is independent of the temperature and adopts the value at the melting point so that the integral is trivial,

Δμ0(T)RT=ΔHmRTm2(TTm),
(10)

where ΔHm=Hwat0Hice0 (determined at Tm) is the melting enthalpy. This equation can be solved in conjunction with Eq. (7) for different values of the water molar fraction xwat, obtaining the freezing temperature at that composition if the solution were ideal,

ΔT=RTm2ΔHmln1xionsRTm2MwΔHmνm.
(11)

In the last step, the logarithm was replaced by a Taylor expansion, valid for low concentrations [ln(1 − xions) ≈ −xions and xionsnions/nwat when xions → 0]. ν is the number of ions in which the salt dissociates. We have also used that m=nsaltnwatMw, nions = nsaltν, and Mw = 0.018 015 kg mol−1 is the molar mass of water. Note that Eq. (11) is the expression mentioned in the Introduction and is often used to estimate the freezing point depression. As we will see later, Eq. (11) provides good estimates of the experimental freezing point depression up to 1 mol kg−1 for 1:1 electrolytes and up to 0.35 mol kg−1 for 1:2 or 2:1 electrolytes.

Alternatively, the integrand of Eq. (9) can be more rigorously evaluated if the dependence of the integrand of Eq. (9) is known. For that aim, we performed NpT simulations of ice and water at different temperatures along the p = 1 bar line. These data were then fitted to a second-degree polynomial. By plugging the fit in Eq. (9), we obtain that for the TIP4P/2005 model of water:

Δμ0(ΔT)RT=9.06×103ΔT+1.33×106ΔT2+4.28×107ΔT3.
(12)

Note that the linear term of this expansion is given by ΔHm/RTm2. Thus, another estimate of the ideal freezing point depression can be obtained by solving Eq. (12) together with Eq. (7) for different values of the water molar fraction xwat.

Finally, the activity of water can be obtained by substituting Eq. (12) in Eq. (6),

lnawat=Δμ0(ΔT)RT.
(13)

The activity of water in an electrolyte solution when in equilibrium with ice depends on the thermodynamic conditions (i.e., p and T) but not on the type of electrolyte in the solution.4,85 The logarithm of the water activity can be expanded as a polynomial function of the freezing point depression,86 whose first coefficient is ΔHm/RTm2.

As at constant temperature and pressure, the chemical potentials of water and salt are related by the Gibbs–Duhem equation,

nwatdμwat+nsaltdμsalt=0,dlnasalt=1mMwdlnawat,
(14)

once the activity of water is known, the activity of the salt can also be obtained86 by integration from very low salt concentrations. Note, however, that for performing this integration, one needs to have accurate measurements of the freezing point depression at very low salt concentrations, something that is not possible using our direct coexistence simulations.

The ice–solution equilibrium curves for NaCl, KCl, LiCl, MgCl2, and Li2SO4 as a function of the salt concentration obtained in this work are given in Fig. 5 and in Table II. Before discussing these, we will show that our results are not significantly affected either by the use of the usual LRC for the LJ potential or by finite size effects. First, to assess the effect of LRC, we repeated the calculations for the NaCl solution starting from a 4 mol kg−1 concentration at T = 236 K using a longer cut-off of 12 Å with LRC. The evolution of the molality of the solution with time is shown in Fig. 6. It can be seen that the equilibrium molality obtained with the 10 Å cut-off [4.07(0.08) mol kg−1] is in agreement with that obtained with the 12 Å cut-off [4.03(0.08) mol kg−1] within the statistical uncertainty of the simulations (see Fig. 5 and Table II). Consistent with this, in previous work, we found that the coexistence between ice and a NaCl solution obtained with a 9 Å cut-off and LRC is the same (within the statistical uncertainty), as that obtained with a cutoff of 13 Å and switching off the LJ LRC.90 Second, to make sure that our results are not significantly affected by finite size effects, we repeated the simulation by doubling the size of the ice slab and the solution along the z-axis for the NaCl system at T = 236 K. In this system, an ice slab with 4096 molecules and 3.6 × 3.1 × 11.7 nm3 dimensions was put in contact with a solution containing 6660 water molecules so that the system dimensions are roughly 3.6 × 3.1 × 30.5 nm3. The results obtained for both system sizes are the same within the statistical uncertainty of the simulations. In particular, the equilibrium concentration was found to be 4.07(0.08) mol kg−1 in the small system and 3.92(0.08) mol kg−1 in the large system. In Table II, mi is the initial concentration of the solution, T is the temperature, m is the average molality at equilibrium, tr is the time the systems are simulated for, and ta is the time accounted to average the concentrations. The uncertainty in the average molality was estimated from the standard deviation of the values used to calculate m.

FIG. 5.

Freezing point depression of ice in contact with (a) NaCl, (b) KCl, (c) LiCl, (d) MgCl2, and (e) Li2SO4 aqueous solutions as a function of salt concentration, as obtained from simulations using the Madrid-2019 model. Experimental data from the literature87–89 (see references therein for the original sources of the experimental measurements) have been fitted to ΔT(m) = −ν 1.856 mbm1.5cm2, and this fit is shown for comparison (blue line). For molalities below 0.1 (1:1) or 0.02 (1:2 and 2:1), the red lines (simulation results) are below the blue lines (experiments) in agreement with the expected limiting behavior at very low concentrations.

FIG. 5.

Freezing point depression of ice in contact with (a) NaCl, (b) KCl, (c) LiCl, (d) MgCl2, and (e) Li2SO4 aqueous solutions as a function of salt concentration, as obtained from simulations using the Madrid-2019 model. Experimental data from the literature87–89 (see references therein for the original sources of the experimental measurements) have been fitted to ΔT(m) = −ν 1.856 mbm1.5cm2, and this fit is shown for comparison (blue line). For molalities below 0.1 (1:1) or 0.02 (1:2 and 2:1), the red lines (simulation results) are below the blue lines (experiments) in agreement with the expected limiting behavior at very low concentrations.

Close modal
TABLE II.

Details of the simulations to estimate the freezing point depression at 1 bar. To assess the effect of changing the model parameters between the Madrid-2019 and the Madrid-2017 model for NaCl, the freezing point was estimated at T = 236 K using both models (the asterisk indicates that the Madrid-2017 model was used). We checked the possible existence of finite size effects by running a simulation with a bigger system (indicated with a b super-index) and the cut-off effect by running a simulation with a cut-off distance of 12 Å (indicated with a R super-index).

Saltmi (mol kg−1)T (K)ΔT (K)m (mol kg−1)tr (ns)ta (ns)
NaCl 247 −3 1.19(0.05) 1750 550 
243 −7 2.13(0.08) 1500 700 
236 −14 4.07(0.08) 1200 500 
4* 236 −14 4.02(0.06) 1400 500 
4b 236 −14 3.92(0.08) 1250 550 
4R 236 −14 4.03(0.08) 1250 550 
224 −26 6.24(0.05) 1100 400 
KCl 247 −3 1.02(0.05) 1250 350 
243 −7 2.20(0.05) 1100 400 
236 −14 3.98(0.05) 900 300 
LiCl 220 −30 4.48(0.04) 2100 400 
225 −25 4.06(0.05) 1200 400 
236 −14 3.30(0.05) 4300 450 
200 −35 6.20(0.08) 2400 800 
MgCl2 232.5 −17.5 2.15(0.05) 1050 750 
215 −35 3.07(0.02) 1150 400 
Li2SO4 233 −17 3.26(0.05) 1250 450 
Saltmi (mol kg−1)T (K)ΔT (K)m (mol kg−1)tr (ns)ta (ns)
NaCl 247 −3 1.19(0.05) 1750 550 
243 −7 2.13(0.08) 1500 700 
236 −14 4.07(0.08) 1200 500 
4* 236 −14 4.02(0.06) 1400 500 
4b 236 −14 3.92(0.08) 1250 550 
4R 236 −14 4.03(0.08) 1250 550 
224 −26 6.24(0.05) 1100 400 
KCl 247 −3 1.02(0.05) 1250 350 
243 −7 2.20(0.05) 1100 400 
236 −14 3.98(0.05) 900 300 
LiCl 220 −30 4.48(0.04) 2100 400 
225 −25 4.06(0.05) 1200 400 
236 −14 3.30(0.05) 4300 450 
200 −35 6.20(0.08) 2400 800 
MgCl2 232.5 −17.5 2.15(0.05) 1050 750 
215 −35 3.07(0.02) 1150 400 
Li2SO4 233 −17 3.26(0.05) 1250 450 
FIG. 6.

Evolution of the molality with time in the NaCl system at T = 236 K in the three setups of the simulation, namely, using the system size adopted in this work (3.6 × 3.1 × 11.7 nm3) using a 10 and a 12 Å cut-off (in both cases using the usual LRC to the LJ contribution), and in a system twice larger in the dimension parallel to the interface (3.6 × 3.1 × 30.5 nm3) using a 10 Å cut-off with the usual LRC to the LJ contribution. The equilibrium values are represented with a solid line.

FIG. 6.

Evolution of the molality with time in the NaCl system at T = 236 K in the three setups of the simulation, namely, using the system size adopted in this work (3.6 × 3.1 × 11.7 nm3) using a 10 and a 12 Å cut-off (in both cases using the usual LRC to the LJ contribution), and in a system twice larger in the dimension parallel to the interface (3.6 × 3.1 × 30.5 nm3) using a 10 Å cut-off with the usual LRC to the LJ contribution. The equilibrium values are represented with a solid line.

Close modal

Let us now look at how the simulated data compares with the available experimental measurements in the literature. To aid that comparison, the simulation freezing points were fitted to a polynomial function of the form

ΔT(m)=ν1.998mbm1.5cm2,
(15)

where m is the molality of the given salt expressed in mol kg−1 and the first coefficient was fixed to the theoretical value RTm2Mw/ΔHm, evaluated using the values for the melting temperature Tm and the melting enthalpy ΔHm of the TIP4P/2005 model.49 The fitted parameters are provided in Table III. Note that for Li2SO4, the fit for the simulated data was not included because the freezing point was calculated for only one salt concentration. As can be seen in Fig. 5, the agreement between the freezing point depression obtained with the Madrid-2019 model and the experimental data87–89 is, in general, in very good agreement for all the salts and at all the concentrations considered (the experimental data are given in the supplementary material). Looking into the results in more detail, it can be concluded that for NaCl, MgCl2, and specially for Li2SO4, the freezing point depression is slightly underestimated with respect to the experimental data over the whole salt concentration range. For KCl and LiCl, the agreement with experiments is rather good for all the salt concentrations.

TABLE III.

Parameters of the quadratic polynomial fit of the freezing point depression of this work’s simulation results, ΔT=TmsolutionTm, as a function of the salt concentration m. The function used for the fit is given by ΔT(m) = −ν 1.998mbm1.5cm2, where m is the molality of the given salt expressed in mol kg−1 and the first coefficient is assigned the theoretical value RTm2Mw/ΔHm=1.998, obtained using the values of Tm and ΔHm for the TIP4P/2005 model, and ν is the number of ions per mol of salt.

Saltbc
NaCl −1.639 0.683 
KCl −1.565 0.665 
LiCl −2.315 1.605 
MgCl2 −6.859 5.676 
Saltbc
NaCl −1.639 0.683 
KCl −1.565 0.665 
LiCl −2.315 1.605 
MgCl2 −6.859 5.676 

As mentioned before, there are not many simulation studies of the freezing point depression that can be used to compare the performance of the Madrid-2019 model with other model potentials commonly used in electrolyte aqueous solutions. We are only aware of three previous estimates of the freezing point depression of water for NaCl solutions, also obtained by using the direct coexistence method, which have been included in Fig. 7 for comparison. The first estimation was performed by Kim and Yethiraj15 who described water using the TIP5P model, and the Åqvist92 and OPLS force fields for Na+ and Cl, respectively. Note that this model uses integer charges for the ions. In the work of Soria et al.,17 the descent in the melting point for a salt solution at a 2 mol kg−1 concentration was estimated using the Joung Cheatham model for NaCl (that was parameterized for SCP/E water and uses integer charges for the ions) and the TIP4P/2005 model for water. This model yields a higher point depression than the Madrid-2019 model, overestimating the experimental result. Conde et al.16 performed a more systematic study, covering several salt concentrations up to 4 mol kg−1. In this case, the solution was modeled using also TIP4P/2005 and a model with integer charges for NaCl.93 The agreement with experiments is comparable to that obtained with the Madrid-2019 model. The difference is that the results of Conde et al. slightly overestimate the experimental results, whereas the Madrid-2019 model underestimates them.

FIG. 7.

Freezing point depression of ice in contact with NaCl aqueous solutions as a function of salt concentration obtained with the Madrid-2019 model, compared to the results for other force fields taken from the literature. The Madrid-201716,17 data were obtained employing force-fields that use TIP4P/2005 for water. For this reason, they provide the same results at low concentrations where the freezing point depression is a colligative property. However, they differ at moderate and high concentrations due to the different force field used for NaCl that provokes that the curve of activity vs concentration is not identical. Kim and Yethiraj employed a force-field that uses TIP5P for water, which predicts a lower value of the constant of proportionality of Eq. (11) (see Table V). This explains why this model underestimates the freezing point depression. The melting point of pure water for TIP5P was taken as 272 K.66Tmsolution taken from Ref. 17 is the one obtained by the chemical potential route, 242 K.

FIG. 7.

Freezing point depression of ice in contact with NaCl aqueous solutions as a function of salt concentration obtained with the Madrid-2019 model, compared to the results for other force fields taken from the literature. The Madrid-201716,17 data were obtained employing force-fields that use TIP4P/2005 for water. For this reason, they provide the same results at low concentrations where the freezing point depression is a colligative property. However, they differ at moderate and high concentrations due to the different force field used for NaCl that provokes that the curve of activity vs concentration is not identical. Kim and Yethiraj employed a force-field that uses TIP5P for water, which predicts a lower value of the constant of proportionality of Eq. (11) (see Table V). This explains why this model underestimates the freezing point depression. The melting point of pure water for TIP5P was taken as 272 K.66Tmsolution taken from Ref. 17 is the one obtained by the chemical potential route, 242 K.

Close modal

Since the Madrid-2019 model modifies the NaCl parameters with respect to a previous model fitted only to properties of this salt (the Madrid-2017 model48) so that the same parameters can be used in several saline solutions, it is pertinent to study whether the freezing point depression is affected by this reparameterization. In principle, since the two model potentials are quite similar, it is expected that both will exhibit the same freezing point depression for the NaCl + H2O system. We checked that this is indeed true by performing additional direct coexistence simulations of the ice–NaCl solution interface at T = 236 K, close to the experimental eutectic point. As can be seen in Fig. 7, the value obtained with the model Madrid-2017 coincides within the statistical uncertainty with that provided by the Madrid-2019 model. Thus, we can conclude that the slight modification of the parameters from the Madrid-2017 to the Madrid-2019 model for NaCl does not change the performance on the estimation of the freezing point depression.

The freezing point depression curves for all the studied salts are gathered in Fig. 8(a), together with the experimental data. This figure illustrates that the Madrid-2019 model provides a good overall description of the effect different salts have on the equilibrium ice–solution curves for all the salts and at all concentrations considered. As described in Sec. II, the freezing point depression, if the mixture were ideal, can be easily computed from the chemical potentials of ice and water in the pure system, and these results have been included in Fig. 8(b) for comparison. Note that there are two ideal curves depending on the number of ions present in the salt (two ions per salt molecule in 1:1 salts and three ions per salt molecule in 2:1 and 1:2 salts). These curves were calculated using the experimental water activity,85 from which it is straightforward to calculate Δμ0T)/RT, and solving Eq. (12) for different values of the water molar fraction xwat. As can be seen in Fig. 8(b), the freezing point exhibits small deviations from the ideal behavior (i.e., lower than 0.4 K) for concentrations lower than 1 mol kg−1 for 1:1 electrolytes and for concentrations lower than 0.35 mol kg−1 for 1:2 or 2:1 electrolytes. Above these concentrations, deviations from ideal behavior are clearly visible.

FIG. 8.

Freezing point depression curves for different salty aqueous solutions. (a) Comparison between the experimental data (solid lines) and the simulation results with the Madrid-2019 model (dashed lines). (b) Comparison between the experimental data (solid lines) and the ideal behavior (dashed lines). The ideal behavior was calculated using the experimental water activity,85 from which it is straightforward to calculate Δμ0T)/RT, and solving Eq. (12) for different values of the water molar fraction xwat. The experimental eutectic points are marked with empty circles.

FIG. 8.

Freezing point depression curves for different salty aqueous solutions. (a) Comparison between the experimental data (solid lines) and the simulation results with the Madrid-2019 model (dashed lines). (b) Comparison between the experimental data (solid lines) and the ideal behavior (dashed lines). The ideal behavior was calculated using the experimental water activity,85 from which it is straightforward to calculate Δμ0T)/RT, and solving Eq. (12) for different values of the water molar fraction xwat. The experimental eutectic points are marked with empty circles.

Close modal

One might think that if the freezing point depression is in good agreement with experiments, this means that the force field should also provide accurate estimates of the activity of water as a function of supercooling and as a function of salt concentration, as well as of the activity of salt (or the salt osmotic coefficient). In what follows, we will show that this is not entirely true.

The natural logarithm of the water activity as a function of supercooling, ln awatT), for the TIP4P/2005 model [calculated using Eq. (13)] is given in Table IV and represented in Fig. 9, along with experimental data taken from Ref. 85. As can be seen in Fig. 9(a), the TIP4P/2005 yields results that are in reasonable agreement with experiments, but there is a small deviation as supercooling increases. How can we explain these differences? As given in Ref. 86, the logarithm of the water activity along the freezing curve can be expanded as a polynomial function of ΔT, whose linear term is given by ΔHm/RTm2, where Tm and ΔHm are, respectively, the melting temperature and the melting enthalpy of pure water. The first term of this expansion is shown in Fig. 9(b) with a dotted line. As can be seen, this term alone already gives a fairly accurate estimation of the water activity, especially in the experimental case in which both lines are practically overlying. For TIP4P/2005, this approximation is also fairly good, but it slightly deteriorates as temperature decreases. From these results, we can conclude that by looking at the value of ΔHm/RTm2, one can get a first indication on whether the water model provides a reasonable estimate of the water activity. The value of this quantity for several popular water models is given in Table V. It seems clear from these data that TIP4P/2005 is the model that gives a better estimation of the water activity, followed by TIP4P/Ice, being TIP5P the one with larger deviations from experiments.

TABLE IV.

Simulation results of the water activities (awat) and activity coefficients (γwat) for the different salts along the freezing depression curve. Δμ0(Tmsolution) is the chemical potential difference between ice and pure water at the freezing point of the considered solution.

Saltm (mol kg−1)Tmsolution (K)Δμ0RTmsolutionawatγwat
NaCl 244.06 −0.054 0.948 0.982 
240.26 −0.088 0.916 0.966 
235.91 −0.126 0.881 0.945 
225.54 −0.215 0.807 0.894 
KCl 243.72 −0.057 0.945 0.979 
240.01 −0.090 0.914 0.963 
235.91 −0.126 0.881 0.945 
LiCl 242.83 −0.065 0.937 0.971 
235.95 −0.126 0.882 0.929 
226.89 −0.204 0.816 0.874 
202.19 −0.385 0.680 0.754 
MgCl2 234.78 −0.136 0.873 0.904 
216.54 −0.286 0.752 0.792 
Saltm (mol kg−1)Tmsolution (K)Δμ0RTmsolutionawatγwat
NaCl 244.06 −0.054 0.948 0.982 
240.26 −0.088 0.916 0.966 
235.91 −0.126 0.881 0.945 
225.54 −0.215 0.807 0.894 
KCl 243.72 −0.057 0.945 0.979 
240.01 −0.090 0.914 0.963 
235.91 −0.126 0.881 0.945 
LiCl 242.83 −0.065 0.937 0.971 
235.95 −0.126 0.882 0.929 
226.89 −0.204 0.816 0.874 
202.19 −0.385 0.680 0.754 
MgCl2 234.78 −0.136 0.873 0.904 
216.54 −0.286 0.752 0.792 
FIG. 9.

(a) Comparison of the natural logarithm of the water activity as a function of the freezing point depression for the TIP4P/2005 water model (red line) and from experiments taken from Ref. 85 (black line). (b) and (c) The natural logarithm of the water activity is compared with (b) an estimate obtained from the first term of the expansion taken from Ref. 86, which is given by (ΔHm/RTm2)ΔT, and (c) with that obtained using the two terms of the expansion given in Ref. 91. In both figures, solid lines are the exact values and the dotted lines are the estimates from the expansions.

FIG. 9.

(a) Comparison of the natural logarithm of the water activity as a function of the freezing point depression for the TIP4P/2005 water model (red line) and from experiments taken from Ref. 85 (black line). (b) and (c) The natural logarithm of the water activity is compared with (b) an estimate obtained from the first term of the expansion taken from Ref. 86, which is given by (ΔHm/RTm2)ΔT, and (c) with that obtained using the two terms of the expansion given in Ref. 91. In both figures, solid lines are the exact values and the dotted lines are the estimates from the expansions.

Close modal
TABLE V.

Evaluation of the factor ΔHm/RTm2 for pure water that provides the linear coefficient of a polynomial expansion of ln awat as a function of the freezing point depression ΔT for different water models typically used to study salt solutions, together with the experimental value. Data for TIP4P/2005 were taken from Ref. 49, for TIP4P/ice from Ref. 69, and for TIP5P from Ref. 64. The last coefficient of the table yields the coefficient of the freezing point depression of an ideal solution [see Eq. (11) of the main text].

ModelTm (K)ΔHm (kcal mol−1)ΔHm/RTm2 (×103 K−1)RTm2Mw/ΔHm (K kg mol−1)
TIP4P/2005 250 1.125 9.06 1.998 
TIP4P/Ice 270 1.29 8.80 2.048 
TIP5P 272 1.75 11.90 1.514 
Exp 273.15 1.44 9.71 1.856 
ModelTm (K)ΔHm (kcal mol−1)ΔHm/RTm2 (×103 K−1)RTm2Mw/ΔHm (K kg mol−1)
TIP4P/2005 250 1.125 9.06 1.998 
TIP4P/Ice 270 1.29 8.80 2.048 
TIP5P 272 1.75 11.90 1.514 
Exp 273.15 1.44 9.71 1.856 

However, coming back to the differences between the experimental activity and that predicted by the TIP4P/2005 model, we have just seen that ΔHm/RTm2 explains much, but not all the discrepancy between experiments and the TIP4P/2005 model, especially as temperature decreases. To explain these differences, we resort to the expansion of the chemical potential difference between ice and water,91 

Δμ0=ΔHmΔT+ΔCp,m(ΔT)2/2Tm,
(16)

where ΔCp,m is the difference between the heat capacity at constant pressure of pure water and ice at the melting point of pure water. Equation (16) can be obtained from the chemical potential difference between ice and water that is obtained when the heat capacity differenceΔCp is assumed to be constant and equal to its value at the melting temperature ΔCp,m,94 

Δμ0=ΔHmTTm1+ΔCp,mTTmT1lnTmT.
(17)

This equation is exact (within the ΔCp = ΔCp,m approximation). Substituting T by Tm + ΔT in Eq. (17) and doing a Taylor expansion with respect to ΔT, one arrives at Eq. (16). Thus, Eq. (16) is approximate, whereas Eq. (17) is exact (within the ΔCp = ΔCp,m approximation).

From Eq. (13), it follows that Δμ0 = RT ln(awat). The heat capacities of ice and water at melting of TIP4P/2005 were taken from Ref. 95, whose values are Cp,ice(T = 250 K) = 13.890 cal mol−1 K−1 and Cp,wat(T = 250 K) = 24.058 cal mol−1 K−1. For the experiments, we used Cp,ice(T = 273 K) = 9.015 cal mol−1 K−196 and Cp,wat(T = 273 K) = 18.013 cal mol−1 K−1.97 ln(awat) estimated from this route is compared with the exact value in Fig. 9(c) and now both the experimental activity and that predicted with TIP4P/2005 are practically coincident with that obtained from Eq. (16) (using, respectively, the experimental and simulated enthalpy of fusion and the heat capacity change at melting). Thus, although the value of the term ΔHm/RTm2 is enough to qualitatively understand the differences between experiments and simulations, the change of heat capacity at melting must also be taken into account to fully understand the discrepancies. Note that ΔCp,m adopts fairly similar values in experiments (ΔCp,m = 8.998 cal mol−1 K−1) and in the TIP4P/2005 model (ΔCp,m = 10.168 cal mol−1 K−1), but the slope of the heat capacity of ice is not well reproduced by classical simulations, as this property is significantly affected by nuclear quantum effects that must be explicitly incorporated in the simulations.95 

Another condition that the force-field should fulfill to reproduce the freezing point depression is that it should describe how the activity of water changes with the concentration of salt. This last dependence is captured by the osmotic pressure, i.e., the pressure at which a solution and the solvent separated by a semipermeable membrane are in equilibrium. Related with the osmotic pressure, one can also measure the osmotic coefficient that can be calculated using the Lewis–Randall definition,

ϕ(T,p,m)=lnawat(T,p,m)νmMwat,
(18)

where ν is the number of ions that result from the dissociation of a salt molecule. In the literature, it is common to find the variation of the osmotic coefficient with the salt concentration at constant temperature and pressure. However, the present simulations do not provide enough information to construct the whole activity curve as a function of both temperature and salt concentration. However, if we combine the activity of water as a function of supercooling, awatT), and the freezing point depression as a function of the salt concentration, ΔT(m), we can estimate the osmotic coefficient of the salt along the freezing curve. Note that awatT) is a universal curve for all the salts,4 but ΔT(m) is different for each salt.

As can been seen in Fig. 10, the osmotic coefficient shows only a mild dependence with temperature. The experimental data at 273 and 298 K are fairly close to each other.98 Coherently with this observation, the osmotic coefficient along the freezing depression curve, obtained from the Lewis–Randall expression in conjunction with the experimental data for awatT) from Ref. 85 and ΔT(m) from Ref. 89, is in good agreement with the curves at constant temperature from Ref. 98, especially with that at T = 273 K up to 3 mol kg−1. Thus, we can conclude that the osmotic coefficient calculated at the freezing depression curves gives a reasonable estimate of the osmotic coefficient at 273 K.

FIG. 10.

Osmotic coefficient of a NaCl solution. (a) Experimental data at T = 273 K and T = 298 K from two different sources (Exp1 from Ref. 98 and Exp2 from Ref. 99) are shown in green lines, finding an excellent agreement between both sets of data. The estimate obtained from the Lewis–Randall expression using the experimental data [awatT) from Ref. 85 and ΔT(m) from Ref. 89] is designated as Exp (freezing line) JCED 2018. The results provided by the Madrid-2019 model are also shown. (b) The osmotic coefficients for other electrolyte models (at T = 298 K taken from Ref. 7) are compared with the Madrid-2019 model and experimental data. JC stands for the non-polarizable Joung–Cheatham force field tailored to SPC/E water.22 AH/SWM4-DP100 and AH/BK3101 are both aqueous alkali halide force fields tailored to the polarizable models of water SWM4-DP102 and BK3.103 For the Madrid-2019 model, besides the results obtained from the fit, the data calculated from the freezing point depressions given in Table II are also shown.

FIG. 10.

Osmotic coefficient of a NaCl solution. (a) Experimental data at T = 273 K and T = 298 K from two different sources (Exp1 from Ref. 98 and Exp2 from Ref. 99) are shown in green lines, finding an excellent agreement between both sets of data. The estimate obtained from the Lewis–Randall expression using the experimental data [awatT) from Ref. 85 and ΔT(m) from Ref. 89] is designated as Exp (freezing line) JCED 2018. The results provided by the Madrid-2019 model are also shown. (b) The osmotic coefficients for other electrolyte models (at T = 298 K taken from Ref. 7) are compared with the Madrid-2019 model and experimental data. JC stands for the non-polarizable Joung–Cheatham force field tailored to SPC/E water.22 AH/SWM4-DP100 and AH/BK3101 are both aqueous alkali halide force fields tailored to the polarizable models of water SWM4-DP102 and BK3.103 For the Madrid-2019 model, besides the results obtained from the fit, the data calculated from the freezing point depressions given in Table II are also shown.

Close modal

Once we confirmed that the osmotic coefficient only exhibits a small dependence with temperature, we again used the Lewis–Randall expression but now with the data obtained from simulations with the Madrid-2019 model for awatT) and ΔT(m). The results are shown in Fig. 10(b), together with the osmotic coefficient for other force-fields at 298 K taken from Ref. 7. The Madrid-2019 model underestimates the osmotic coefficient up to about a 20%, but the predicted curve runs almost parallel to the experimental one above a 1 mol kg−1 concentration. We have included in Fig. 10(b) the osmotic coefficient evaluated using the larger system. Although its value is slightly closer to the experimental one, the deviation is still large, showing that finite size effects are not responsible for the deviation between experiment and the value of the osmotic coefficient predicted by the Madrid-2019 force field. The osmotic coefficient is related to the variation of the chemical potential of water when salt is added.

The osmotic coefficients of the other force-fields included in Fig. 10 exhibit even larger departures from experiment, failing to reproduce the correct dependence with salt concentration. This includes the popular Joung–Cheetham force-field (that uses SCP/E for water) and two polarizable force-fields, all with integer charges. These models correctly predict the osmotic coefficient at low concentrations but fail to capture the dependence with salt concentration. Thus, it seems that although the Madrid-2019 model does not give the correct absolute value (due to the use of scaled charges), it provides the best dependence with concentration.

In Secs. IV A and IV B, we have seen that the Madrid-2019 model does not give an accurate description of either the water activity or the osmotic coefficient, so how is it possible that this model gives a reasonable description of the freezing point depression? The reason is that there is a cancellation of errors between the water activity and the salt osmotic coefficient. The overestimation of the water activity of the TIP4P/2005 model is compensated by the description of the impact of adding salt to the solution.

Now that we have better understood what the contributions of the water activity and the salt osmotic coefficient to the estimation of the freezing point depression are, it is pertinent to analyze what the desired behavior for a good model of electrolyte solutions should be. As we have just seen, some of the deficiencies of the electrolytes force-fields describing the freezing point depression arise from deficiencies of the considered water model to reproduce the water activity. The changes in the chemical potential of salt and the activity of water are related by the Gibbs–Duhem relation [Eq. (14)]. If a force-field correctly reproduced the water activity for any value of m, that means that it would also reproduce the osmotic coefficient and the activity coefficient of the salt, as the osmotic coefficient and the ln γ are related to each other at constant temperature and pressure.84 Making the approximation that the changes in the osmotic coefficient (and, consequently, in the logarithm of the activity coefficient of the salt, ln γ) with temperature at constant composition are small, we can calculate the target freezing point in the following way. First, for a given salt molality, m, we look for the experimental freezing point depression [ΔTexp in Fig. 9(a)]. Second, for this value of ΔTexp, we look for the experimental value of ln awat [Fig. 9(a)]. Third, we see for which ΔT, the water model gives the experimental ln awat, and this would be the target freezing point depression ΔTtarget for molality m [Fig. 9(a)]. In this way, we are subtracting the inherent error coming from the inability of the water model to give the right values of the water activity. This target curve is approximate, since it was obtained by disregarding the variation of ϕ or ln γ with temperature. However, as can be seen by comparing the experimental osmotic coefficient of NaCl along the freezing curve and at T = 273 K shown in Fig. 10(a), this is a fairly good approximation up to 4 mol kg−1 and reasonable beyond that concentration. The target freezing depression curve for NaCl using the TIP4P/2005 model for water is given in Fig. 11. If we want to have a model for NaCl solutions using TIP4P/2005 that correctly describes the salt osmotic coefficient, we should aim to reproduce this target curve, rather than the true experimental one. In this case, when fitting or testing the model parameters, one should aim to overestimate the freezing point depression following the target curve. Note that this curve is universal for any NaCl solution model that uses TIP4P/2005 to describe water. Similar target curves have been estimated for other salts [KCl in Fig. 11(b) and LiCl in Fig. 11(c)]. Note that the larger departures of the target from the experimental curve are observed for LiCl, as this salt is the one for which the eutectic point appears at higher salt concentrations.

FIG. 11.

Target freezing point depression curves for the NaCl (a), KCl (b) and LiCl (c) solutions for a model that properly reproduces the change of the activity of water as a function of the salt concentration. These curves take into account that the TIP4P/2005 model does not correctly reproduce the activity of water as a function of supercooling.

FIG. 11.

Target freezing point depression curves for the NaCl (a), KCl (b) and LiCl (c) solutions for a model that properly reproduces the change of the activity of water as a function of the salt concentration. These curves take into account that the TIP4P/2005 model does not correctly reproduce the activity of water as a function of supercooling.

Close modal

In this work, we have shown that the Madrid-2019 model provides a faithful description of the experimental freezing point depression of five salt solutions (NaCl, KCl, LiCl, MgCl2, and Li2SO4) at concentrations up to the eutectic point. The best agreement is obtained for KCl and LiCl, whereas the freezing point depression of NaCl, MgCl2, and Li2SO4, is slightly underestimated over the whole concentration range that covers up to the eutectic point. By analyzing the possible causes of the differences between the simulated and experimental data, we discovered that the water model TIP4P/2005 introduces some error in the water activity along the freezing curve. Thus, when fitting or testing the goodness of an electrolyte model, it is wise not to try to reproduce the experimental freezing point, but rather one should target a freezing depression curve obtained by subtracting the error introduced by the use of a water model that does not reproduce exactly the activity of water. Otherwise, the electrolyte model can only provide an accurate estimation of the freezing point depression by a cancellation of errors between the description of the water activity and the salt osmotic coefficient. We have also seen that the Madrid-2019 model underestimated the osmotic coefficient of NaCl over the whole concentration range, but at molalities higher than 1 mol kg−1, the simulated and experimental curves run parallel to each other. This underestimation of the osmotic coefficient compensates to some extent the higher activity of water of the TIP4P/2005 model for a certain supercooling, thus yielding reasonable freezing curves due to a certain cancellation of errors. In the future, it would be of interest to check the performance of other force fields of electrolytes for predicting the freezing point depression. Note that reproducing the freezing point depression without cancellation of errors would require a water model with the right melting temperature and melting enthalpy and a force field of ions able to reproduce the osmotic coefficient up to high concentrations for temperatures well below the melting temperature of pure water. Certainly a force field like that has not been proposed so far. The Madrid-2019 force field is far from being perfect, but it does a reasonable job in describing the experimental results. This adds to number of properties whose description improves with the Madrid-2019 model as compared to other non-polarizable models with integer charges, which include the densities, viscosities, and diffusion coefficients of salt solutions at atmospheric conditions and up to the solubility limit,33 and even the solubility limit of NaCl.48 The trade-off is that these type of potentials are not adequate to describe properties of pure salts8,57 (either in solid or liquid state), in which the charge distribution of the ions is expected to be very different from that in aqueous solution.

We hope the freezing point depression becomes a standard property to be tested when proposing force fields for electrolytes as its determination is relatively straightforward and feasible with current computational resources. In addition, it provides easy and straightforward information about the osmotic coefficient electrolyte solutions. The osmotic coefficient is calculated easily but unfortunately not along an isotherm but along the freezing point depression curve.

See the supplementary material for an ASCII with the experimental freezing point depression data obtained from Refs. 87–89 (see references therein for the original sources of the experimental measurements).

This project was funded by Grant No. PID2019-105898GB-C21 of the Ministry of Science, Innovation and Universities. E.G.N. acknowledges Agencia Estatal de Investigación and Fondo Europeo de Desarrollo Regional (FEDER) (Grant No. PID2020-115722GB-C21). C.P.L. acknowledges Ministerio de Universidades for a predoctoral Formacion Profesorado Universitario (Grant No. FPU18/03326) and also Ayuntamiento de Madrid for a Residencia de Estudiantes grant. The authors acknowledge the computer resources and technical assistance provided by the Red Española de Supercomputación (RES).

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article.

1.
L.
Vrbka
and
P.
Jungwirth
, “
Brine rejection from freezing salt solutions: A molecular dynamics study
,”
Phys. Rev. Lett.
95
,
148501
(
2005
).
2.
V. F.
Petrenko
and
R. W.
Whitworth
,
Physics of Ice
(
Oxford University Press
,
Oxford
,
1999
).
3.
P. G.
Debenedetti
,
Metastable Liquids: Concepts and Principles
(
Princeton University Press
,
1996
).
4.
T.
Koop
,
B.
Luo
,
A.
Tsias
, and
T.
Peter
, “
Water activity as the determinant for homogeneous ice nucleation in aqueous solutions
,”
Nature
406
,
611
(
2000
).
5.
D.
Frenkel
and
A. J. C.
Ladd
, “
New Monte Carlo method to compute the free energy of arbitrary solids. Application to the fcc and hcp phases of hard spheres
,”
J. Chem. Phys.
81
,
3188
3193
(
1984
).
6.
J. L.
Aragones
,
E. G.
Noya
,
C.
Valeriani
, and
C.
Vega
, “
Free energy calculations for molecular solids using GROMACS
,”
J. Chem. Phys.
139
,
034104
(
2013
).
7.
W. R.
Smith
,
F.
Moučka
, and
I.
Nezbeda
, “
Osmotic pressure of aqueous electrolyte solutions via molecular simulations of chemical potentials: Application to NaCl
,”
Fluid Phase Equilib.
407
,
76
(
2016
).
8.
A. Z.
Panagiotopoulos
, “
Simulations of activities, solubilities, transport properties, and nucleation rates for aqueous electrolyte solutions
,”
J. Chem. Phys.
153
,
010903
(
2020
).
9.
A. J. C.
Ladd
and
L. V.
Woodcock
, “
Triple-point coexistence properties of the Lennard-Jones system
,”
Chem. Phys. Lett.
51
,
159155
(
1977
).
10.
C.
Vega
,
E.
Sanz
,
J. L. F.
Abascal
, and
E. G.
Noya
, “
Determination of phase diagrams via computer simulation: Methodology and applications to water, electrolytes and proteins
,”
J. Phys.: Condens. Matter
20
,
153101
(
2008
).
11.
L.
Vrbka
and
P.
Jungwirth
, “
Molecular dynamics of freezing of water and salt solutions
,”
J. Mol. Liq.
134
,
64
(
2007
).
12.
M. A.
Carignano
,
E.
Baskaran
,
P. B.
Shepson
, and
I.
Szleifer
, “
Molecular dynamics simulation of ice growth from supercooled pure water and from salt solution
,”
Ann. Glaciol.
44
,
113
(
2006
).
13.
I.
Tsironi
,
D.
Schlesinger
,
A.
Späh
,
L.
Eriksson
,
M.
Segad
, and
F.
Perakis
, “
Brine rejection and hydrate formation upon freezing of NaCl aqueous solutions
,”
Phys. Chem. Chem. Phys.
22
,
7625
(
2020
).
14.
S.
Luo
,
Y.
Jin
,
R.
Tao
,
H.
Li
,
C.
Li
,
J.
Wang
, and
Z.
Li
, “
Molecular understanding of ion rejection in the freezing of aqueous solutions
,”
Phys. Chem. Chem. Phys.
23
,
13292
13299
(
2021
).
15.
J. S.
Kim
and
A.
Yethiraj
, “
The effect of salt on the melting of ice: A molecular dynamics simulation study
,”
J. Chem. Phys.
129
,
124504
(
2008
).
16.
M. M.
Conde
,
M.
Rovere
, and
P.
Gallo
, “
Molecular dynamics simulations of freezing-point depression of TIP4P/2005 water in solution with NaCl
,”
J. Mol. Liq.
261
,
513
519
(
2018
).
17.
G. D.
Soria
,
J. R.
Espinosa
,
J.
Ramirez
,
C.
Valeriani
,
C.
Vega
, and
E.
Sanz
, “
A simulation study of homogeneous ice nucleation in supercooled salty water
,”
J. Chem. Phys.
148
,
222811
(
2018
).
18.
J. R.
Espinosa
,
C.
Vega
,
C.
Valeriani
, and
E.
Sanz
, “
Seeding approach to crystal nucleation
,”
J. Chem. Phys.
144
,
034501
(
2016
).
19.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
, “
Comparison of simple potential functions for simulating liquid water
,”
J. Chem. Phys.
79
,
926
935
(
1983
).
20.
H. W.
Horn
,
W. C.
Swope
,
J. W.
Pitera
,
J. D.
Madura
,
T. J.
Dick
,
G. L.
Hura
, and
T.
Head-Gordon
, “
Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew
,”
J. Chem. Phys.
120
,
9665
(
2004
).
21.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
, “
The missing term in effective pair potentials
,”
J. Phys. Chem.
91
,
6269
6271
(
1987
).
22.
I. S.
Joung
and
T. E.
Cheatham
, “
Determination of alkali and halide monovalent ion parameters for use in explicit solvated biomolecular simulation
,”
J. Phys. Chem. B
112
,
9020
(
2008
).
23.
I.
Nezbeda
,
F.
Moučka
, and
W. R.
Smith
, “
Recent progress in molecular simulation of aqueous electrolytes: Force fields, chemical potentials and solubility
,”
Mol. Phys.
114
,
1665
1690
(
2016
).
24.
W. R.
Smith
,
I.
Nezbeda
,
J.
Kolafa
, and
F.
Moučka
, “
Recent progress in the molecular simulation of thermodynamic properties of aqueous electrolyte solutions
,”
Fluid Phase Equilib.
466
,
19
30
(
2018
).
25.
Z.
Mester
and
A. Z.
Panagiotopoulos
, “
Temperature-dependent solubilities and mean ionic activity coefficients of alkali halides in water from molecular dynamics simulations
,”
J. Chem. Phys.
143
,
044505
(
2015
).
26.
Z.
Mester
and
A. Z.
Panagiotopoulos
, “
Mean ionic activity coefficients in aqueous NaCl solutions from molecular dynamics simulations
,”
J. Chem. Phys.
142
,
044507
(
2015
).
27.
J. R.
Espinosa
,
J. M.
Young
,
H.
Jiang
,
D.
Gupta
,
C.
Vega
,
E.
Sanz
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
, “
On the calculation of solubilities via direct coexistence simulations: Investigation of NaCl aqueous solutions and Lennard-Jones binary mixtures
,”
J. Chem. Phys.
145
,
154111
(
2016
).
28.
N. F. A.
van der Vegt
,
K.
Haldrup
,
S.
Roke
,
J.
Zheng
,
M.
Lund
, and
H. J.
Bakker
, “
Water-mediated ion pairing: Occurrence and relevance
,”
Chem. Rev.
116
,
7626
7641
(
2016
).
29.
C. J.
Fennell
,
A.
Bizjak
,
V.
Vlachy
, and
K. A.
Dill
, “
Ion pairing in molecular simulations of aqueous alkali halide solutions
,”
J. Phys. Chem. B
113
,
6782
6791
(
2009
).
30.
Y.
Yao
and
Y.
Kanai
, “
Free energy profile of NaCl in water: First principles molecular dynamics with SCAN and wB97X-V exchange-correlation functionals
,”
J. Chem. Theory Comput.
14
,
884
893
(
2018
).
31.
A. K.
Giri
and
E.
Spohr
, “
Cluster formation of NaCl in bulk solutions: Arithmetic vs. geometric combination rules
,”
J. Mol. Liq.
228
,
63
(
2017
).
32.
S.
Yue
and
A. Z.
Panagiotopoulos
, “
Dynamic properties of aqueous electrolyte solutions from non-polarisable, polarisable, and scaled-charge models
,”
Mol. Phys.
117
,
3538
3549
(
2019
).
33.
I. M.
Zeron
,
J. L. F.
Abascal
, and
C.
Vega
, “
A force field of Li+, Na+, K+, Mg2+, Ca2+, Cl, and SO42 in aqueous solution based on the TIP4P/2005 water model and scaled charges for the ions
,”
J. Chem. Phys.
151
,
134504
(
2019
).
34.
J. S.
Kim
,
Z.
Wu
,
A. R.
Morrow
,
A.
Yethiraj
, and
A.
Yethiraj
, “
Self-diffusion and viscosity in electrolyte solutions
,”
J. Phys. Chem. B
116
,
12007
12013
(
2012
).
35.
Z. R.
Kann
and
J. L.
Skinner
, “
A scaled-ionic-charge simulation model that reproduces enhanced and suppressed water diffusion in aqueous salt solutions
,”
J. Chem. Phys.
141
,
104507
(
2014
).
36.
I. V.
Leontyev
and
A. A.
Stuchebrukhov
, “
Electronic continuum model for molecular dynamics simulations
,”
J. Chem. Phys.
130
,
085102
(
2009
).
37.
I. V.
Leontyev
and
A. A.
Stuchebrukhov
, “
Electronic polarizability and the effective pair potentials of water
,”
J. Chem. Theory Comput.
6
,
3153
3161
(
2010
).
38.
I.
Leontyev
and
A.
Stuchebrukhov
, “
Accounting for electronic polarization in non-polarizable force fields
,”
Phys. Chem. Chem. Phys.
13
,
2613
2626
(
2011
).
39.
E.
Pluhařová
,
P. E.
Mason
, and
P.
Jungwirth
, “
Ion pairing in aqueous lithium salt solutions with monovalent and divalent counter-anions
,”
J. Phys. Chem. A
117
,
11766
11773
(
2013
).
40.
M.
Kohagen
,
P. E.
Mason
, and
P.
Jungwirth
, “
Accurate description of calcium solvation in concentrated aqueous solutions
,”
J. Phys. Chem. B
118
,
7902
7909
(
2014
).
41.
M.
Kohagen
,
P. E.
Mason
, and
P.
Jungwirth
, “
Accounting for electronic polarization effects in aqueous sodium chloride via molecular dynamics aided by neutron scattering
,”
J. Phys. Chem. B
120
,
1454
1460
(
2015
).
42.
E.
Duboué-Dijon
,
P. E.
Mason
,
H. E.
Fischer
, and
P.
Jungwirth
, “
Hydration and ion pairing in aqueous Mg2+ and Zn2+ solutions: Force-field description aided by neutron scattering experiments and ab initio molecular dynamics simulations
,”
J. Phys. Chem. B
122
,
3296
3306
(
2017
).
43.
T.
Martinek
,
E.
Duboué-Dijon
,
Š.
Timr
,
P. E.
Mason
,
K.
Baxová
,
H. E.
Fischer
,
B.
Schmidt
,
E.
Pluhařová
, and
P.
Jungwirth
, “
Calcium ions in aqueous solutions: Accurate force field description aided by ab initio molecular dynamics and neutron scattering
,”
J. Chem. Phys.
148
,
222813
(
2018
).
44.
M.
Předota
and
D.
Biriukov
, “
Electronic continuum correction without scaled charges
,”
J. Mol. Liq.
314
,
113571
(
2020
).
45.
J.
Li
and
F.
Wang
, “
Pairwise-additive force fields for selected aqueous monovalent ions from adaptive force matching
,”
J. Chem. Phys.
143
,
194505
(
2015
).
46.
E. E.
Bruce
and
N. F. A.
van der Vegt
, “
Does an electronic continuum correction improve effective short-range ion-ion interactions in aqueous solution?
,”
J. Chem. Phys.
148
,
222816
(
2018
).
47.
R.
Fuentes-Azcatl
and
M. C.
Barbosa
, “
Sodium chloride, NaCl/ɛ: New force field
,”
J. Phys. Chem. B
120
,
2460
2470
(
2016
).
48.
A. L.
Benavides
,
M. A.
Portillo
,
V. C.
Chamorro
,
J. R.
Espinosa
,
J. L. F.
Abascal
, and
C.
Vega
, “
A potential model for sodium chloride solutions based on the TIP4P/2005 water model
,”
J. Chem. Phys.
147
,
104501
(
2017
).
49.
J. L. F.
Abascal
and
C.
Vega
, “
A general purpose model for the condensed phases of water: TIP4P/2005
,”
J. Chem. Phys.
123
,
234505
(
2005
).
50.
C.
Vega
, “
Water: One molecule, two surfaces, one mistake
,”
Mol. Phys.
113
,
1145
(
2015
).
51.
M.
Jorge
and
L.
Lue
, “
The dielectric constant: Reconciling simulation and experiment
,”
J. Chem. Phys.
150
,
084108
(
2019
).
52.
I. M.
Zeron
,
M. A.
Gonzalez
,
E.
Errani
,
C.
Vega
, and
J. L. F.
Abascal
, “
‘In silico’ seawater
,”
J. Chem. Theory Comput.
17
,
1715
1725
(
2021
).
53.
S.
Blazquez
,
M. M.
Conde
,
J. L. F.
Abascal
, and
C.
Vega
, “
The Madrid-2019 force field for electrolytes in water using TIP4P/2005 and scaled charges: Extension to the ions F, Br, I, Rb+, and Cs+
,”
J. Chem. Phys.
156
,
044505
(
2022
).
54.
R. C.
Weast
,
M. J.
Astle
, and
W. H.
Beyer
,
CRC Handbook of Chemistry and Physics
(
CRC Press
,
Boca Raton, FL
,
1986
).
55.
T.
Yagasaki
,
M.
Matsumoto
, and
H.
Tanaka
, “
Lennard-Jones parameters determined to reproduce the solubility of NaCl and KCl in SPC/E, TIP3P, and TIP4P/2005 water
,”
J. Chem. Theory Comput.
16
,
2460
(
2020
).
56.
S.
Blazquez
,
I. M.
Zeron
,
M. M.
Conde
,
J. L. F.
Abascal
, and
C.
Vega
, “
Scaled charges at work: Salting out and interfacial tension of methane with electrolyte solutions from computer simulations
,”
Fluid Phase Equilib.
513
,
112548
(
2020
).
57.
D.
Kussainova
,
A.
Mondal
,
J. M.
Young
,
S.
Yue
, and
A. Z.
Panagiotopoulos
, “
Molecular simulation of liquid–vapor coexistence for NaCl: Full-charge vs scaled-charge interaction models
,”
J. Chem. Phys.
153
,
024501
(
2020
).
58.
Y.
Ding
,
A. A.
Hassanali
, and
M.
Parrinello
, “
Anomalous water diffusion in salt solutions
,”
Proc. Natl. Acad. Sci. U. S. A.
111
,
3310
3315
(
2014
).
59.
Y.
Yao
,
M. L.
Berkowitz
, and
Y.
Kanai
, “
Communication: Modeling of concentration dependent water diffusivity in ionic solutions: Role of intermolecular charge transfer
,”
J. Chem. Phys.
143
,
241101
(
2015
).
60.
B.
Doherty
,
X.
Zhong
,
S.
Gathiaka
,
B.
Li
, and
O.
Acevedo
, “
Revisiting OPLS force field parameters for ionic liquid simulations
,”
J. Chem. Theory Comput.
13
,
6131
6145
(
2017
).
61.
M.
Soniat
and
S. W.
Rick
, “
The effects of charge transfer on the aqueous solvation of ions
,”
J. Chem. Phys.
137
,
044511
(
2012
).
62.
M.
Soniat
,
G.
Pool
,
L.
Franklin
, and
S. W.
Rick
, “
Ion association in aqueous solution
,”
Fluid Phase Equilib.
407
,
31
38
(
2016
).
63.
J.
Kolafa
, “
Solubility of NaCl in water and its melting point by molecular dynamics in the slab geometry and a new BK3-compatible force field
,”
J. Chem. Phys.
145
,
204509
(
2016
).
64.
C.
Vega
,
J. L. F.
Abascal
,
M. M.
Conde
, and
J. L.
Aragones
, “
What ice can teach us about water interactions: A critical comparison of the performance of different water models
,”
Faraday Discuss.
141
,
251
276
(
2009
).
65.
C.
Vega
,
M.
Martin-Conde
, and
A.
Patrykiejew
, “
Absence of superheating for ice Ih with a free surface: A new method of determining the melting point of different water models
,”
Mol. Phys.
104
,
3583
3592
(
2006
).
66.
R.
García Fernández
,
J. L. F.
Abascal
, and
C.
Vega
, “
The melting point of ice Ih for common water models calculated from direct coexistence of the solid-liquid interface
,”
J. Chem. Phys.
124
,
144506
(
2006
).
67.
M. M.
Conde
,
M.
Rovere
, and
P.
Gallo
, “
High precision determination of the melting points of water TIP4P/2005 and water TIP4P/Ice models by the direct coexistence technique
,”
J. Chem. Phys.
147
,
244506
(
2017
).
68.
M. W.
Mahoney
and
W. L.
Jorgensen
, “
A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions
,”
J. Chem. Phys.
112
,
8910
8922
(
2000
).
69.
J. L. F.
Abascal
,
E.
Sanz
,
R.
García Fernández
, and
C.
Vega
, “
A potential model for the study of ices and amorphous water: TIP4P/Ice
,”
J. Chem. Phys.
122
,
234511
(
2005
).
70.
A. J. C.
Ladd
and
L. V.
Woodcock
, “
Interfacial and co-existence properties of the Lennard-Jones system at the triple point
,”
Mol. Phys.
36
,
611
(
1978
).
71.
E. G.
Noya
,
C.
Vega
, and
E.
de Miguel
, “
Determination of the melting point of hard spheres by direct coexistence simulations
,”
J. Chem. Phys.
128
,
154507
(
2008
).
72.
J. D.
Bernal
and
R. H.
Fowler
, “
A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions
,”
J. Chem. Phys.
1
,
515
(
1933
).
73.
V.
Buch
,
P.
Sandler
, and
J.
Sadlej
, “
Simulations of H2O solid, liquid, and clusters, with an emphasis on ferroelectric ordering transition in hexagonal ice
,”
J. Phys. Chem. B
102
,
8641
(
1998
).
74.
D.
Rozmanov
and
P. G.
Kusalik
, “
Temperature dependence of crystal growth of hexagonal ice (Ih)
,”
Phys. Chem. Chem. Phys.
13
,
15501
(
2011
).
75.
M. M.
Conde
,
C.
Vega
, and
A.
Patrykiejew
, “
The thickness of a liquid layer on the free surface of ice as obtained from computer simulations
,”
J. Chem. Phys.
129
,
014702
(
2008
).
76.
B.
Hess
,
C.
Kutzner
,
D.
van der Spoel
, and
E.
Lindahl
, “
GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation
,”
J. Chem. Theory Comput.
4
,
435
447
(
2008
).
77.
J. R.
Espinosa
,
E.
Sanz
,
C.
Valeriani
, and
C.
Vega
, “
On fluid-solid direct coexistence simulations: The pseudo-hard sphere model
,”
J. Chem. Phys.
139
,
144502
(
2013
).
78.
M.
Parrinello
and
A.
Rahman
, “
Polymorphic transitions in single crystals: A new molecular dynamics method
,”
J. Appl. Phys.
52
,
7182
7190
(
1981
).
79.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
, “
Nosé–Hoover chains: The canonical ensemble via continuous dynamics
,”
J. Chem. Phys.
97
,
2635
2643
(
1992
).
80.
D. R.
Wheeler
and
J.
Newman
, “
A less expensive Ewald lattice sum
,”
Chem. Phys. Lett.
366
,
537
543
(
2002
).
81.
B.
Hess
,
H.
Bekker
,
H. J. C.
Berendsen
, and
J. G. E. M.
Fraaije
, “
LINCS: A linear constraint solver for molecular simulations
,”
J. Comput. Chem.
18
,
1463
1472
(
1997
).
82.
B.
Hess
, “
P-LINCS: A parallel linear constraint solver for molecular simulation
,”
J. Chem. Theory Comput.
4
,
116
122
(
2008
).
83.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J. C.
Berendsen
, “
Numerical integration of the Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes
,”
J. Comput. Phys.
23
,
327
341
(
1977
).
84.
I. M.
Klotz
and
R. M.
Rosenberg
,
Chemical Thermodynamics. Basic Concepts and Methods
(
John Willey and Sons
,
NJ
,
2008
).
85.
H.
Sippola
and
P.
Taskinen
, “
Activity of supercooled water on the ice curve and other thermodynamic properties of liquid water up to the boiling point at standard pressure
,”
J. Chem. Eng. Data
63
,
2986
(
2018
).
86.
I. M.
Klotz
and
R. M.
Rosenberg
,
Termodinámica Química: Teoría y Métodos Básicos
(
Editorial AC
,
Madrid
,
1977
).
87.
A. N.
Campbell
, “
The system Li2SO4-H2O
,”
J. Am. Chem. Soc.
65
,
2268
2271
(
1943
).
88.
D.
Li
,
D.
Zeng
,
H.
Han
,
L.
Guo
,
X.
Yin
, and
Y.
Yao
, “
Phase diagrams and thermochemical modeling of salt lake brine systems. I. LiCl + H2O system
,”
Calphad
51
,
1
12
(
2015
).
89.
D.
Li
,
D.
Zeng
,
X.
Yin
,
H.
Han
,
L.
Guo
, and
Y.
Yao
, “
Phase diagrams and thermochemical modeling of salt lake brine systems. II. NaCl + H2O, KCl + H2O, MgCl2 + H2O and CaCl2 + H2O systems
,”
Calphad
53
,
78
89
(
2016
).
90.
V.
Bianco
,
M. M.
Conde
,
C. P.
Lamas
,
E. G.
Noya
, and
E.
Sanz
, “
Phase diagram of the NaCl-water system from computer simulations
,”
J. Chem. Phys.
156
,
064505
(
2022
).
91.
C.
Wang
,
J.
Wu
,
H.
Wang
, and
Z.
Zhang
, “
Classical nucleation theory of ice nucleation: Second-order corrections to thermodynamic parameters
,”
J. Chem. Phys.
154
,
234503
(
2021
).
92.
J.
Åqvist
, “
Ion-water interaction potentials derived from free energy perturbation simulations
,”
J. Phys. Chem.
94
,
8021
(
1990
).
93.
M. A.
Portillo
,
J. L. F.
Abascal
, and
C.
Vega
, “
Simulación de disoluciones acuosas por ordenador
,” M.Sc. thesis,
Universidad Complutense de Madrid
,
2016
.
94.
M. A. R.
Martins
,
S. P.
Pinho
, and
J. A. P.
Coutinho
, “
Insights into the nature of eutectic and deep eutectic mixtures
,”
J. Solution Chem.
48
,
962
982
(
2019
).
95.
C.
Vega
,
M. M.
Conde
,
C.
McBride
,
J. L. F.
Abascal
,
E. G.
Noya
,
R.
Ramirez
, and
L. M.
Sesé
, “
Heat capacity of water: A signature of nuclear quantum effects
,”
J. Chem. Phys.
132
,
046101
(
2010
).
96.
R.
Feistel
and
W.
Wagner
, “
A new equation of state for H2O ice Ih
,”
J. Phys. Chem. Ref. Data
35
,
1021
(
2006
).
97.
C. A.
Angell
,
D. R.
MacFarlane
, and
M.
Oguni
, “
The Kauzmann paradox, metastable liquids, and ideal glasses: A summary
,”
Ann. N. Y. Acad. Sci.
484
,
241
(
1986
).
98.
L. J.
Partanen
,
J. I.
Partanen
, and
P. O.
Minkkinen
, “
Traceable values for activity and osmotic coefficients in aqueous sodium chloride solutions at temperatures from 273.15 to 373.15 K up to the saturated solutions
,”
J. Chem. Eng. Data
65
,
5226
5239
(
2020
).
99.
D. G.
Archer
, “
Thermodynamic properties of the NaCl + H2O system. II. Thermodynamic properties of NaCl(aq), NaCl·2H2(cr), and phase equilibria
,”
J. Phys. Chem. Ref. Data
21
,
793
(
1992
).
100.
G.
Lamoureux
and
B.
Roux
, “
Absolute hydration free energy scale for alkali and halide ions established from simulations with a polarizable force field
,”
J. Phys. Chem. B
110
,
3308
3322
(
2006
).
101.
P. T.
Kiss
and
A.
Baranyai
, “
A new polarizable force field for alkali and halide ions
,”
J. Chem. Phys.
141
,
114501
(
2014
).
102.
G.
Lamoureux
,
A. D.
MacKerell
, Jr.
, and
B.
Roux
, “
A simple polarizable model of water based on classical Drude oscillators
,”
J. Chem. Phys.
119
,
5185
5197
(
2003
).
103.
P. T.
Kiss
and
A.
Baranyai
, “
A systematic development of a polarizable potential of water
,”
J. Chem. Phys.
138
,
204507
(
2013
).

Supplementary Material