Molecular dynamics simulation of the surface tension of aqueous sodium chloride: from dilute to highly supersaturated solutions and molten salt

Sodium chloride (NaCl) is one of the key components of atmospheric aerosols. The surface tension of aqueous NaCl solution (σNaCl,sol) and its concentration dependence are essential to determine the equilibrium water vapor pressure of aqueous NaCl droplets. Supersaturated NaCl solution droplets are observed in laboratory experiments and under atmospheric conditions, but the experimental data for σNaCl,sol are mostly limited up to subsaturated solutions. In this study, the surface tension of aqueous NaCl is investigated by molecular dynamics (MD) simulations and the pressure tensor method from dilute to highly supersaturated solutions. We show that the linear approximation of concentration dependence of σNaCl,sol at molality scale can be extended to the supersaturated NaCl solution until a molality of ∼ 10.7 molkg−1 (i.e., solute mass fraction (xNaCl) of ∼ 0.39). Energetic analyses show that this monotonic increase in surface tension is driven by the increase in excess surface enthalpy (1H ) as the solution becomes concentrated. After that, the simulated σNaCl,sol remains almost unchanged until xNaCl of ∼ 0.47 (near the concentration upon efflorescence). The existence of the “inflection point” at xNaCl of ∼ 0.39 and the stable surface tension of xNaCl between ∼ 0.39 and ∼ 0.47 can be attributed to the nearly unchanged excess surface entropy term (T ·1S) and the excess surface enthalpy term (1H ). After a “second inflection point” at xNaCl of ∼ 0.47, the simulated σNaCl,sol gradually regains the growing momentum with a tendency to approach the surface tension of molten NaCl (∼ 175.58 mNm−1 at 298.15 K, MD simulation-based extrapolation). This fast increase in σNaCl,sol at xNaCl > 0.47 is a process driven by excess surface enthalpy and excess surface entropy. Our results reveal different regimes of concentration dependence of the surface tension of aqueous NaCl at 298.15 K: a water-dominated regime (xNaCl from 0 to ∼ 0.39), a transition regime (xNaCl from ∼ 0.39 to ∼ 0.47) and a molten NaCl-dominated regime (xNaCl from ∼ 0.47 to 1).


Introduction
Sodium chloride (NaCl) is one of the most important components of atmospheric aerosol particles (Finlayson-Pitts, 2003;Lewis and Schwartz, 2004). The aqueous NaCl solution droplet could participate in many atmospheric processes, such as phase transition, cloud activation, ice crystallization, long-range transport and chemical aging (Martin, 2000;Finlayson-Pitts, 2003;Ghorai et al., 2014;Wagner et al., 2015;Chen et al., 2016). To better understand these processes, the concentration-dependent surface tension of aqueous NaCl solution (σ , ) is essential to determine the equilibrium between NaCl solution droplet and water vapor (Jarvis and Scheiman, 1968;Dutcher et 5 al., 2010).
Because of the energy barrier of crystallization during dehydration and size-effects at nanoscale 10 (Martin, 2000;Biskos et al., 2006;Cheng et al., 2015), supersaturated aqueous NaCl solution droplets can exist under atmospheric conditions. However, direct measurements of surface tension of supersaturated droplets are challenging due to technical difficulties (Harkins and Brown, 1919;Vargaftik et al., 1983;Richardson and Snyder, 1994;Kumar, 2001). Only recently, Bzdek et al. (2016) overcame this limitation with an optical tweezer method and extend the concentration range to ~8 mol 15 kg -1 , where the near linear relationship still holds (Bzdek et al., 2016).
It is a matter of debate to which extent the approximation of a near linear dependence of surface tensions on molality can still be used for NaCl droplets. Cheng et al. (2015) used the Differential Köhler Analyses (DKA) method to retrieve the surface tension of NaCl aqueous droplets, and revealed a large deviation from the near linear increase at molality of ~10 mol kg -1 . In literature, such deviation 20 in concentrated solution has also been found for other compounds, such as HNO 3 (Weissenborn and Pugh, 1996) and it is believed to be typically true for most highly soluble electrolytes (Dutcher et al., 2010). The reason for such deviation remains unclear.
A few surface tension models have been developed for highly concentrated electrolyte solutions, e.g., Li and Lu (2001), Li et al. (1999), Levin and Flores-Mena (2001). Li and Lu (2001) developed a 25 model based on the Gibbs dividing surface concept, where the adsorption and desorption rate constants, saturated surface excess, stoichiometric coefficient of ions and mean ionic activity coefficient are needed. For NaCl aqueous solution, this model is suitable for solution with concentration up to ~5.5 mol kg -1 . Li et al. (1999) uses Debye-Huckel parameter, osmotic coefficient and a proportionality constant from the fitting of measured values to calculate the surface tension, which covers the 30 concentration until saturation point of bulk NaCl aqueous solutions. The remaining models are mostly only suitable for the dilute electrolyte solutions, such as the one proposed by Levin and Flores-Mena (2001). In their valid concentration range, these surface tension models produce linear or near linear concentration dependence of σ , that agrees well with currently available observations. One surface tension model that is able to predict σ , in the whole concentration range from 35 infinitely dilute ( = 0) to highly supersaturated solution to molten salts ( = 1) was proposed by Dutcher et al. (2010), which has been adopted into the widely used Extended Aerosol Inorganics Model (E-AIM) (Wexler and Clegg, 2002). This model is based on the concept: ions are solvated by the water at low salt concentrations, meaning that water-structures like hydration shells are formed around the ions; while at very high salt concentration the water is considered as "solute" hat is solvated by the ions instead, meaning a salt-structure is formed around the water molecules. Accordingly, for a diluted solution, the surface tension of water dominates and the surface tension of the solution equals the surface tension of water adjusted by a term that is proportional to the solute concentration.

MD simulation
MD simulations were carried out with the GROMACS 5.1 package (Abraham et al., 2015). The Na + ions, Clions and water molecules were added into a cubic box ( = 5 nm) to imitate the NaCl solution. The concentrations of simulated solutions are summarized in Table 1. To simulate the surface tension of supersaturated NaCl aqueous solution, we make use of the time window in the MD 25 simulations before the crystallization starts in the system. The highest we can reach is ~0.64 (the corresponding concentration is ~30.39 mol kg -1 ), below which the simulated surface tensions in three independent runs stably converge after 50 to 100 ns ( Fig. 1). For more concentrated solutions, stable convergence cannot be reached, as for example large fluctuations are shown in Fig. 1d  The procedure (Fig. 2) of simulation we followed is: (1) systems were firstly energetically 35 minimized by the steepest-descent method (Stillinger and Weber, 1985) ( 2013) because the system that we were dealing with is much more concentrated. 1 fs time step was adopted and conformations for analysis were saved every 2 ps. Both electrostatic interactions and van der Waals interactions were calculated using the particle mesh Ewald (PME) algorithm, which has been proven to be a good choice for accurate calculation of long-range interactions (Essmann et al., 1995;Fischer et al., 2015). To test the reproducibility, all the systems were simulated 3 times, and the 10 respective statistical error bars were provided.
In our simulation, the Joung-Cheatham (JC) force field for NaCl (Joung and Cheatham III, 2009) with SPC/E water model (Berendsen et al., 1987)

Calculation of Surface Tension
Based on results from MD simulations, the surface tension was calculated by using the mechanical definition of the atomic pressure (Alejandre et al., 1995): where σ can represent the surface tension of molten NaCl (σ ), NaCl solution (σ , ) or pure water (σ ), is the length of the simulation cell in the longest direction (along z-axis) and (a=x, y, z) denotes the diagonal component of the pressure tensor. The 〈… 〉 refers to the time average.
The factor 0.5 outside the squared brackets takes into account the two interfaces in the system. Only the 25 stable values were taken as our calculated surface tension.

Energy analysis
The excessive surface enthalpy denotes the additional enthalpy in the system due to the creation of surfaces. It can be calculated as the difference of enthalpy between solutions with and without surfaces 30 (Bahadur et al., 2007), where _ is the total enthalpy of simulated systems with surfaces and is the total enthalpy of simulated systems without surfaces. As the kinetic energy is the same for systems with or without surfaces and the difference of pV can be ignored, ∆H can be presented as where _ and are the potential energy of the system with and without surfaces.
Then the surface tension can be determined by the excessive surface free energy per unit area as in Eq.
(4) (Davidchack and Laird, 2003): where ∆ is the increased part of free energy due to the creation of surfaces, is the total area of the surface we created, so = 2 × and is the area of each created surface. ∆ is the excessive surface entropy. We then can retrieve ∆ by using the data of enthalpy and surface tension: 5 ∆H and T · ∆S per unit area ( ∆ and ·∆ ) are obtained as the enthalpic and entropic part of contributions to the net surface tension, which will be used to explain the change of surface tension along with the mass fraction of NaCl ( ).

10
3.1 Water-dominated regime ( < ~0.36) In Fig. 3a, the calculated surface tension of NaCl aqueous solution (σ , ) are compared with experimentally determined values (Jarvis and Scheiman, 1968;Johansson and Eriksson, 1974;Aveyard and Saleem, 1976;Weissenborn and Pugh, 1995;Matubayasi et al., 2001;Pegram andRecord, 2006, 2007;Morris et al., 2015;Bzdek et al., 2016)  ) from models and experiments converge nicely (Fig. 3b). The MD simulation is able to reproduce the increment in the growth of surface tension from pure water due to the addition of solute NaCl though the predicted absolute value of σ , is systematically underestimated, which may mainly be attributed to the discrepancy between observed σ and the modeled ones from the 30 SPE/C water model.
By performing MD simulations in the supersaturated concentration range, we found that this linear relationship still holds beyond the saturation point until of ~ 0.36 (Fig. 4). As mentioned above, the laboratory experiments with elevated NaCl aqueous droplet and the optical tweezer method show that the linear relationship between σ , and NaCl concentration (molality scale) can be extended 35 to ~ 8 mol kg -1 (Fig. 3) (Bzdek et al., 2016), corresponding to of ~ 0.33 (Fig. 4), which is consistent with our simulations.
3.2 Transition regime ( from ~0.36 to ~0.47) Surface tensions of single inorganic electrolyte aqueous solution were often found to be linear functions of concentration (at the molality scale) over moderate concentration range (Talbot, 1987).
However, these simple relationships will not hold when the solutions become more concentrated. As shown in Fig. 4 (Cheng et al., 2015), where a large deviation of surface tension from the monotonic linear increase. Note that beyond of ~0.47, the simulated surface tension increases again (Fig. 4). This "second inflection point", right at the concentration upon efflorescence, may imply potential correlation with crystallization processes. Beyond the "second inflection point" ( > 0.47), the simulated σ , gradually regains a growing momentum (Fig. 4). Unfortunately, due to the large fluctuation in the surface tension simulation (Fig. 1), we are not able to extend our surface tension calculation in this way beyond 15 of ~0.64. However, according to Dutcher et al. (2010), it is expected that the surface tension of the solution would ultimately approach the surface tension of the hypothetical molten solute (i.e., = 1) at the same temperature. This hypothesis has been found to be consistent with the DKA retrieval for a highly concentrated ammonium sulfate aqueous solution with molality of ~380 mol kg -1 (Cheng et al., 2015). We thus also try to constrain the growth of σ   (Frenken and Van der Veen, 1985), which sometimes is observed when the temperature is raised towards the triple point. In our case, the enrichment zone of NaCl (which is about  Fig.7) would be a precursor effect to the (metastable) NaCl-rich bulk solution.
As shown in Fig. 6