Freezing point depression of salt aqueous solutions using the Madrid-2019 model

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 best^1,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,³ 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).⁴ 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 methods^5,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 solution^1,11 (paying special attention to brine rejection from the ice), and this has been followed by other authors.^12–14 In 2008, Kim and Yethiraj¹⁵ 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 depression¹⁶ 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.¹⁷ 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.,¹⁷ kinetic aspects of the nucleation of ice in salty water were also studied, and nucleation rates were estimated using the seeding technique.¹⁸ 

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,¹⁹ TIP4P-Ew,²⁰ or SPC/E.²¹ 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.²² 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.⁸ 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.³⁴ 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 Skinner³⁵ 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 Stuchebrukhov^36–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-workers^39–43 and also implemented by several groups.^44–47 Following this line of research in 2017, we proposed a force field for NaCl in water⁴⁸ (with water described by the TIP4P/2005 model⁴⁹) 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,³⁵ 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 model⁴⁸ is estimated to be 5.7 mol kg⁻¹ at room conditions, in reasonable agreement with the experimental solubility, 6.1 mol kg⁻¹.⁵⁴ 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.⁵⁵ 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 energies⁵⁶ 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 liquids⁶⁰) 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.⁵⁹ 

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 Yethiraj¹⁵ 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.²⁷ to determine the solubility of a salt and by Kolafa⁶³ 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,⁶⁴ and for SPC/E, it is 215 K.⁶⁵ 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 model²⁰ with a melting point of about 244 K⁶⁶ and an even better choice would be the TIP4P/2005 with a melting temperature of 250 K.⁶⁷ 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. TIP5P⁶⁸ and/or TIP4P/Ice⁶⁹ 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 stable⁶⁴ 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,

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

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.⁷¹ 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.⁷¹ 

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.²⁷ 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.¹⁶ 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.³³ This force field describes water and the sulfate ions as rigid and non-polarizable. Water is modeled with the TIP4P/2005 potential⁴⁹ 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 I_h 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 I_h with proton disorder and almost zero dipole moment satisfying the Bernal–Fowler⁷² rules was generated using the algorithm proposed by Buch et al.⁷³ Initial solution concentrations between 1 and 6 mol kg⁻¹ were considered to cover temperatures within the range T = 200–247 K. The corresponding numbers of ions (Li⁺, Na⁺, K⁺, Mg²⁺, 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 L_x ≈ 3.6 nm, L_y ≈ 3.1 nm, and L_z ≈ 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 nm³. Note that the width of the ice slab is large enough to avoid finite size effects on the direct coexistence method.²⁷ 

The system is then allowed to evolve. The simulations are performed using the molecular dynamics GROMACS package⁷⁶ 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,⁷⁸ with a relaxation time of 2 ps and a compressibility of 10⁻⁶. The temperature was set using the Nosé–Hoover thermostat⁷⁹ 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 summations⁸⁰ 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 SHAKE⁸³ 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.

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 = L_z/100 and measuring the water or salt density in each of them,

where w_i(z + Δz) is the mass of the water molecules (w_wat) or of the ions (w_ions) whose positions along the z-coordinate are between z − Δz/2 and z + Δz/2 and ⟨L_x⟩ and ⟨L_y⟩ 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 MgCl₂ 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

where M_salt is the molar mass of the ionic compound. As density profiles are measured in kg/m³, the molar mass M_salt must be expressed in kg mol⁻¹ 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 MgCl₂ 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⁻¹ 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.

The profiles of the molarities of the ions (mol dm⁻³) 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.²⁷ 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.

When comparing with experiments, the melting temperature will be expressed as $ΔT=Tmsolution−Tm$⁠, where T_m 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,⁶⁷ 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 (a_wat),

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 x_wat 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

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

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

where n_wat is the number of moles of water and n_ions,+ and n_ions,− 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

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, T_m, 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,

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

In the last step, the logarithm was replaced by a Taylor expansion, valid for low concentrations [ln(1 − x_ions) ≈ −x_ions and x_ions ≈ n_ions/n_wat when x_ions → 0]. ν is the number of ions in which the salt dissociates. We have also used that $m=nsaltnwatMw$⁠, n_ions = n_saltν, and M_w = 0.018 015 kg mol⁻¹ 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⁻¹ for 1:1 electrolytes and up to 0.35 mol kg⁻¹ 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:

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 x_wat.

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

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,⁸⁶ 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,

once the activity of water is known, the activity of the salt can also be obtained⁸⁶ 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, MgCl₂, and Li₂SO₄ 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⁻¹ 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⁻¹] is in agreement with that obtained with the 12 Å cut-off [4.03(0.08) mol kg⁻¹] 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.⁹⁰ 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 nm³ 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 nm³. 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⁻¹ in the small system and 3.92(0.08) mol kg⁻¹ in the large system. In Table II, m_i is the initial concentration of the solution, T is the temperature, m is the average molality at equilibrium, t_r is the time the systems are simulated for, and t_a 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.

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

where m is the molality of the given salt expressed in mol kg⁻¹ and the first coefficient was fixed to the theoretical value $RTm2Mw/ΔHm$⁠, evaluated using the values for the melting temperature T_m and the melting enthalpy ΔH_m of the TIP4P/2005 model.⁴⁹ The fitted parameters are provided in Table III. Note that for Li₂SO₄, 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 data^87–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, MgCl₂, and specially for Li₂SO₄, 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.

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 Yethiraj¹⁵ who described water using the TIP5P model, and the Åqvist⁹² 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.,¹⁷ the descent in the melting point for a salt solution at a 2 mol kg⁻¹ 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.¹⁶ performed a more systematic study, covering several salt concentrations up to 4 mol kg⁻¹. In this case, the solution was modeled using also TIP4P/2005 and a model with integer charges for NaCl.⁹³ 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.

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 model⁴⁸) 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 + H₂O 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,⁸⁵ from which it is straightforward to calculate Δμ⁰(ΔT)/RT, and solving Eq. (12) for different values of the water molar fraction x_wat. 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⁻¹ for 1:1 electrolytes and for concentrations lower than 0.35 mol kg⁻¹ for 1:2 or 2:1 electrolytes. Above these concentrations, deviations from ideal behavior are clearly visible.

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 a_wat(ΔT), 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 T_m and ΔH_m 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.

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,⁹¹ 

where ΔC_p,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ΔC_p is assumed to be constant and equal to its value at the melting temperature ΔC_p,m,⁹⁴ 

This equation is exact (within the ΔC_p = ΔC_p,m approximation). Substituting T by T_m + Δ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 ΔC_p = ΔC_p,m approximation).

From Eq. (13), it follows that Δμ⁰ = RT ln(a_wat). The heat capacities of ice and water at melting of TIP4P/2005 were taken from Ref. 95, whose values are C_p,ice(T = 250 K) = 13.890 cal mol⁻¹ K⁻¹ and C_p,wat(T = 250 K) = 24.058 cal mol⁻¹ K⁻¹. For the experiments, we used C_p,ice(T = 273 K) = 9.015 cal mol⁻¹ K⁻¹⁹⁶ and C_p,wat(T = 273 K) = 18.013 cal mol⁻¹ K⁻¹.⁹⁷ ln(a_wat) 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 ΔC_p,m adopts fairly similar values in experiments (ΔC_p,m = 8.998 cal mol⁻¹ K⁻¹) and in the TIP4P/2005 model (ΔC_p,m = 10.168 cal mol⁻¹ K⁻¹), 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.⁹⁵ 

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,

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, a_wat(ΔT), 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 a_wat(ΔT) is a universal curve for all the salts,⁴ 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.⁹⁸ 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 a_wat(ΔT) 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⁻¹. 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.

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 a_wat(ΔT) 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⁻¹ 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.⁸⁴ 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 [ΔT_exp in Fig. 9(a)]. Second, for this value of ΔT_exp, we look for the experimental value of ln a_wat [Fig. 9(a)]. Third, we see for which ΔT, the water model gives the experimental ln a_wat, and this would be the target freezing point depression ΔT_target 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⁻¹ 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.

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, MgCl₂, and Li₂SO₄) at concentrations up to the eutectic point. The best agreement is obtained for KCl and LiCl, whereas the freezing point depression of NaCl, MgCl₂, and Li₂SO₄, 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⁻¹, 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,³³ and even the solubility limit of NaCl.⁴⁸ The trade-off is that these type of potentials are not adequate to describe properties of pure salts^8,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.