1 Introduction

In many industrial processes concentrated salt solutions are of absolute necessity. One such process is the production of sodium chlorate that requires concentrated salt mixtures up to 6 mol·L−1 [1]. Concentrated solutions are also found in nature such as naturally occurring brines. In order to understand the reactions occurring in these concentrated electrolytes, basic thermodynamic data such as activity and osmotic coefficients are needed. Such data can be acquired experimentally or through calculations. The Debye–Hückel (DH) theory of electrolytes provides simple analytical equations, which can be used to calculate the activity and osmotic coefficients [2]. However, the activity coefficients can only be calculated up to about 0.01 mol·L−1 due to the neglect of ion size, short range interactions and ion–ion correlations. Over the years there have been numerous efforts to extend the applicability range of the DH theory by including the above mentioned effects empirically, i.e., by use of fitting parameters. The best system of equations available at present is the Pitzer model by which activity and osmotic coefficients of single as well as mixed salts at different temperatures can be modelled up to ionic strengths as high as 6 mol·L−1 [3]. It is not surprising that the Pitzer model is widely used in industrial applications because of the complexity of the electrolyte systems. Although the Pitzer model is very useful it requires quite a number of fitting parameters, therefore physical transparency of the fitted parameters can be lost.

Monte Carlo (MC) simulation is an alternative method to the theoretical approach and is comparatively less approximate. To date, the MC method has not been fully exploited in the field of electrolytes since it requires a substantial amount of computer resources. However, with the rapid development in computer speed and capacity, one can foresee that simulations may become a method of choice when dealing with complex electrolytes [4]. In recent studies, the feasibility of the MC method to predict the activity and osmotic coefficients of single salts at 25 °C [5] as well as salt mixtures such as sea water in the temperature range of 0 to 25 °C has been shown [6, 7]. These and other studies [8, 9] showed that by adjusting only the ion size, the experimentally determined activity and osmotic coefficients can be modelled to concentrations much higher than the DH theory. Moreover, it was shown that the fitted ionic radii obtained from single salts can be used to simulate the properties of salt mixtures [4, 5].

In recent years molecular dynamics (MD) simulations have also been used to calculate the osmotic, activity coefficients [10,11,12] and the solubility of salt solutions [13, 14]. Although progress has been made in calculating the osmotic and activity coefficients by MD simulations that are comparable with the experimental data, the results of MD simulations are highly dependent on the force field used in the simulations and on the system’s size [10, 15]. Especially, simulations of bulk properties at the molecular level, such as the concentration dependent dielectric constant of a salt, is challenging as shown by Hess et al. [10].

Another advantage of MC simulations is that single ion activity coefficients are calculated and mean ionic activity coefficients can be obtained by taking the average. At present there is no reliable theory by which single ion activities can be calculated. A long time ago Bates [16] provided a system of equations by which mean ionic activity coefficients could be split into single ion activity coefficients by using hydration numbers. Later on, Ferse [17] developed a system of equations by which single ion activities can be calculated from mean values. The most recent work is of Fraenkel [18], who developed the smaller ion shell (SIS) model based on the DH theory to calculate single ion activity coefficients. It is interesting to note that the single ion, i.e. cation or anion, activity coefficients calculated by the above-mentioned theories do not agree with each other, probably due to the varying approximations used in these theories. On the experimental side there have been a number of reports in recent literature where single ion activities measured by ion selective electrodes have been reported [19,20,21,22,23]. However, the measurements of single ion activities by ion sensitive electrodes has been criticized in the literature [24, 25]. The criticism was based on the one of pioneering comments made by Guggenheim [26] stating that “single ion activities cannot be measured” and the debate whether they can be measured or not by ion selective electrodes is still going on [27,28,29]. One major obstacle in measuring single ion activity coefficients by ion selective electrodes is the interdependence of activity coefficients and liquid junction potential as recently discussed by Ferse et al. [30]. Considering the different views on the calculations and measurements of single ion activities it is very tempting and interesting to compare the single ion activity coefficients calculated by MC simulations with corresponding values obtained by ion selective electrodes as well as those calculated by different theories. Since MC simulations are less approximate than other theories, we may be able to shed some light on this old controversy.

The article is structured as follows. Feasibility of the MC method for the simulation of mean ionic activity and osmotic coefficients of industrial relevant salts by using the experimentally determined dielectric constant of salt solutions at various concentrations is shown in Sects. 3.1 and 3.2, and for the activity coefficients of acids and salt solutions up to 70 °C by using the temperature dependent dielectric constant of water is shown in Sect. 3.3. Sections 3.4 and 3.5 describe the feasibility of the MC method for calculating the water activities in concentrated salt solutions; pH and ionic activities in mixed salt solutions is also considered. Finally, in Sect. 3.6 results regarding the comparison of single ion activity coefficients obtained by MC with experiment and other theoretically predicted values are discussed.

2 Methods

2.1 Theoretical

In this study MC simulations were performed with the unrestricted primitive model (UPM) of electrolytes. In the UPM, cations and anions are considered as charged hard spheres of different radii. The solvent is considered as a continuum having a specified dielectric constant which for water at 25 °C is 78.36. Let us consider a cation as a sphere with radius ri and an anion of radius rj. The interaction between two charged spheres, U(rij), is described by Eq. 1 (Fig. 1):

Fig. 1
figure 1

A schematic representation of the electrolyte model used in MC simulations

Full size image
$$\begin{aligned} u\left( {r_{ij} } \right) = \frac{{z_{i} z_{j} e^{2} }}{{4\pi \varepsilon_{0} \varepsilon_{s} r_{ij} }}\quad r_{ij} > r_{i} + r_{j} \hfill \\ u\left( {r_{ij} } \right) = \infty \quad r_{ij} \le r_{i} + r_{j} \hfill \\ \end{aligned}$$
(1)

Here zi and zj are the charge of ions (e.g., + 1 and – 1 for Na+ and Cl−), e is the electron charge, ε0 is permittivity of vacuum and εs is the permittivity of the solvent (in this case water). This model depicts that when the distance between the two charged ions (rij) is larger than the sum of their radii then the electrostatic interaction will be given by the screened columbic potential, whereas when (rij) is less than the sum of the individual radii the interaction between the charged particles will become infinite. This is because the spheres are considered as hard spheres that are non-overlapping.

The above mentioned electrostatic model was implemented in the MC simulation method. The MC simulations were performed by the standard Metropolis algorithm. A canonical ensemble was constructed by using a cubic box with periodic boundary conditions. In order to attain a required concentration of electrolyte in the cubic box the side length was varied for a fixed number of total particles. To obtain activity coefficients, a modified Widom particle insertion method was used [31]. The Widom method states that a non-overlapping particle α inserted at a random position r will have the activity coefficient γ given by Eq. 2:

$$\ln \gamma = - \ln \left\langle {\exp \left[ { - \beta \Delta U_{\alpha } \left( r \right)} \right]} \right\rangle_{0}$$
(2)

The exponential term enclosed in angled brackets is the ensemble average of the energy change, ΔUα, of adding a particle and β = 1/kBT. This method provides a direct calculation of the chemical potential. However, the original Widom method becomes less accurate when dealing with ionic systems of finite sizes, since the addition of a charged particle will violate electroneutrality in the cell. This effect can be reduced considerably by using very large systems, but such simulations require enormous computation times. Svensson and Woodward [32], proposed a simple charge rescaling method to re-establish electroneutrality in the computational cell. This method has shown good results for symmetric as well as for asymmetric electrolytes [33] and was used in the present work.

In the MC simulations reported here, large systems, i.e., 1000 ions for 1:1, and 2000–2310 ions for mixtures of electrolytes were used. System size dependence of the calculated activity and osmotic coefficients was evaluated by performing the MC simulations with varying configurations. The number of MC configurations chosen to generate the thermodynamic data was based on the criterion that the variation in the calculated values was minimum, i.e., values up to three decimal places did not vary with the number of MC configurations. The results reported in this study are based on long simulations, i.e., 95 million MC configurations were performed to generate the thermodynamic data. Activity and osmotic coefficients of salts calculated by the MC simulations were fitted to the experimental data. A best fit between experiment and MC was obtained by adjusting the cation radius while keeping the anion radius equal to its crystallographic value.

2.2 Conversion of Experimental and Theoretical Data

Experimentally, the thermodynamic properties of solutions such as vapor pressure are usually determined at molal concentrations (mole of solute per kilogram of solvent), whereas theoretical calculations are performed at molar concentrations. In order to inter-compare the theoretical data with experimental data, the molal scale needs to be converted to the molar scale or vice versa. In this work, conversions from the molal to molar scale was done by using the density of solution according to Eqs. 3 and 4 as described in detail in our previous papers [5, 8].

$$\ln \rho_{s} = \ln \rho_{0} + \frac{AmW}{mW + 1000}$$
(3)
$$c = \frac{{m\rho_{s} }}{1 + 0.001mW}$$
(4)

where ρ0 is the density of pure water at 25 °C, ρs is the density of the solution at 25 °C (g·cm−3), A is a fitting parameter, W is the molecular weight of the electrolyte, m is the molal concentration and c is the molar concentration. Note that to calculate the density of solution by using Eq. 3, the value of an empirical constant A is required. The values of constant A used for calculating the densities of LiClO3 and NaClO3 are 0.64436 and 0.6930, respectively [34]. The conversions of activity and osmotic coefficients from the molar to molal scale were also done by using the density of the solution. A detailed description and equations used for such conversion are given in our previous paper [9].

3 Results and Discussion

The fitting procedure adopted in this study, i.e., the best fitting of calculated activity or osmotic coefficients to the experimental data by adjusting the cation radius while keeping the anion radius fixed at crystallographic value, has its merits in experimental evidence. It has been shown by a number of experimental techniques [35]. especially dielectric spectroscopy [36], that large anions such as Cl−, I−, ClO −4 and ClO −3 are weekly hydrated and can be considered as unhydrated ions. On the other hand, cations have small size and concentrated charge density, which results in strong hydration. Thus, the fitting approach adopted in this study is an attempt to take into account the hydration effects by using the fitted cation radii. This implies that if the fitted cation radius is larger than the crystallographic radius then the difference is due to hydration. We keep in mind that this is a simple approach based on only one fitting parameter, whereas in reality the situation is much more complex. For example, as the salt concentration increases it is possible that ionic hydration layers overlap [37] and in many salts ion pairs are formed [38]. At present such effects are not included in the electrolyte model used for MC simulations but can indirectly be observed by deviations between experiment and simulations as described below. Generally for a good fit between the experimental and theoretical data, deviations for the activity and osmotic coefficients are less than 0.03 and 0.04, respectively.

3.1 Ion Pairing or Ion Dissociation

In salt solutions where ion pairing is present the mean ionic activity coefficients are low lying, compared with salts where the ions are dissociated. The curves describing the mean ion activity and osmotic coefficients at increasing salt concentration bend upward at lower salt concentrations for dissociated salts due to strong ion–ion repulsions. These effects can be captured to some extent in the MC simulations because the ions are modelled as charged hard spheres of different sizes. By comparing the calculated mean ionic activity and osmotic coefficients with the experimental data, a specific concentration above which ion pairs dominate in solution can be localized. Above that concentration the calculated mean ionic activity and osmotic coefficients will have higher values than the experimental values because in the MC simulations ions are modelled as non-overlapping charged hard spheres. In our recent work the feasibility of the MC method to accurately model the activity and osmotic coefficients for a large number of salts was shown [5]. In the present study we have extended the approach to chlorate salts.

Calculated activity and osmotic coefficients of NaClO3 and LiClO3 are compared with the experimental data [34] in Fig. 2. In the case of NaClO3 the best fit was restricted to 1 mol·L−1 whereas in the case of LiClO3 the experimental data were fitted to slightly higher concentration and the quality of the fit was better. However, the most striking differences between these salts were the fitted cation radii. The best-fit cation radius of Li+ is 0.21 nm whereas that of Na+ is 0.102 nm, with the radius of ClO −3 being kept constant at its crystallographic value, i.e., 0.20 nm. The Li+ radius found in LiClO3 is the same that was found for LiCl, i.e., 0.21 nm [5]. On the other hand, the Na+ radius in NaClO3 is equal to the crystallographic value, 0.102 nm) while in NaCl it is 0.168 nm, which is much larger than its crystallographic value. If we calculate the hydration number from the best fitted radius of Li+, as was done in our previous study [5], it is between 3 and 4 while for Na+ it is 0. The lower activity coefficients of NaClO3 and smaller fitted radius of Na+ compared to NaCl indicate that in NaClO3 solutions the Na+ ions are associated with the chlorate ion (ClO −3 ) [39]. On the other hand, the small, highly charged, Li+ ion retains its hydration structure and does not strongly associate with the chlorate ion [36, 40]. This contrasts with recent experimental studies using large angle x-ray scattering and EXAFS [40], where chlorate and perchlorate ions were found to be associated with three water molecules per oxygen, i.e., 9 water molecules for chlorate and 12 for perchlorate. In all cases, the sodium ion was found to be six-coordinated with water. However, it is well known that results obtained with different techniques are not easily compared [36]. The hydration number observed depends on the strength of the ion–water interaction in comparison to the water–water interactions and, even though water molecules are structured around for example chlorate, the resulting evolution of the activity coefficient and osmotic pressure with increasing concentration is best described by a naked anion and varying coordinated water molecules depending on the size of the cation.

Fig. 2
figure 2

a Mean ionic activity and b osmotic coefficients calculated by MC simulations are compared with experimentally determined data. Exp denotes the experimental data [34] and concentration (c) is given as square root of mol·L−1. MC simulations were performed at 25 °C, the dielectric constant of water was 78.36 and ionic radii of Li+, Na+ and ClO −3 were fixed at 0.21, 0.102 and 0.20 nm, respectively, for MC simulations

Full size image

Recent MD simulations [41] in molten NaClO3 show that both Li+ and Na+ make bi-dentate complexes, in agreement with ab initio calculations in the absence of water [42] and IR spectra [43] of matrix isolated Li+–ClO −3 and Na+–ClO −3 ion pairs.

In the chlorate production process, highly concentrated electrolyte mixtures of the major components Na+, ClO −3 , Cl−, ClO −4 , SO 2−4 and chromate ions, in that order, are used. The major cationic impurities are Ca2+ and K+. One can envision that Ca2+ in NaClO3 solutions will behave like the Li+ ion because the Ca2+ ion is strongly hydrated and has a similar hydration structure as the Li+ ion, i.e., two hydration shells [35]. On the other hand, weakly hydrated K+ will be more strongly associated to the chlorate ion than the Na+ ion. Knowledge about the ion–ion interactions and possible ion-pair formation is important for chlorate production, in particular for estimating the solubility of the salt [39] and to optimize the process.

3.2 Activity Coefficients by Using a Concentration Dependent Dielectric Constant

Concentration dependent dielectric constants have been measured for a number of salt solutions [44]. Generally, the dielectric constant of water decreases as the concentration of ions increases. However, the extent to which it changes depends upon the salt type. The concentration dependent dielectric constant has been used as a fitting parameter in theories such as mean spherical approximation (MSA) to calculate the activity coefficients of different salt solutions [45]. Good agreement between the calculated and measured activity coefficients was obtained up to concentrations much higher than for the case of fixed dielectric constant. In a recent study we have shown that the applicability range is extended up to 3 mol·L−1 for KI and KBr solutions by using experimentally determined values in the Corrected Debye–Hückel theory (CDH) and MC simulations instead of a fixed value [9]. In the present study additional salts have been investigated by using the concentration dependent dielectric constant in MC simulations and the results are shown in Fig. 3.

Fig. 3
figure 3

Mean ionic activity coefficients, calculated by MC simulations using concentration dependent dielectric constants of LiCl and NaClO4 from [44] and of NaCl from [46], are compared with the experimentally determined data [47]. Exp denotes experimental data and the concentration (c) is given as square root of mol·L−1. MC simulations were performed at 25 °C and the ionic radii of Li+, Na+ in LiCl, NaCl and NaClO4 were fixed at 0.248, 0.202 and 0.143 nm, respectively. The ionic radii of Cl− and ClO4− were fixed at 0.181 and 0.240 nm, respectively; the anionic radii are the crystallographic radii

Full size image

Note that experimental data for the concentration dependent dielectric constant of LiCl and NaClO4 solutions were taken from Barthel’s book [44] whereas for NaCl they are from Buchner [46]. The experimental data of NaCl was fitted by the polynomial as given in Eq. 5:

$$y(c) = y(0) - \delta c + \beta c^{{\frac{3}{2}}} : \, y = \varepsilon ,\tau$$
(5)

The values of the coefficients y(0), δ, and β are given in Ref. [46] for calculating the concentration dependent dielectric constant ε(c) and relaxation time τ. Similarly, the equations used to calculate the concentration dependent dielectric constant of LiCl and NaClO4 can be found in our previous article [9] or in Barthel’s book [44].

The reason was that the activity coefficients of NaCl calculated by using the Barthel data are much lower than experimental values at concentrations higher than 1 mol·L−1. This inconsistency of NaCl data from Barthel [44] has also been pointed out in a recent study [48]. In the case of LiCl the inclusion of a concentration dependent dielectric constant resulted in good agreement between the calculated and measured activity coefficients up to 3.5 mol·L−1 and for NaCl and NaClO4 good fits up to 2 mol·L−1 were obtained. Compared to the case in which the dielectric constant was fixed [5], this is a considerable increase in concentration range. The best fitted ionic radii of Li+ and Na+ in LiCl, NaCl and NaClO4 are 0.248, 0.202 and 0.143 nm, respectively. These fitted radii are larger than the crystallographic radii indicating strongly hydrated cations in these salt solutions, with on average 5.7, 2.7 and 0.7 water molecules, respectively. These results are consistent with the effective hydration numbers obtained by dielectric spectroscopy measurements [36] and the number of transported water molecules obtained from conductivity measurements combined with modelling [49]. For example the effective hydration numbers for Li+ and Na+ in chloride salts are 8–12, and 4.5, respectively, whereas the number of water molecules transported by Li+ and Na+ are 6 and 3.6, respectively. In both these studies it is pointed out that for anions, except for very small OH− and F− ions, the effective hydration number is zero and number of transported water molecules is close to zero. Thus, the MC results presented in this study (Fig. 3) indicate that a good agreement with the experimental data can be obtained to high concentrations by using the experimentally determined dielectric constants of salt solutions with only using one adjustable parameter, i.e., the cation radius.

3.3 Mixed Electrolytes

Calculation of ionic activity coefficients in mixed electrolyte solutions is much more complicated than for single electrolytes. In theories a number of cross terms are needed to accurately calculate the ion activity coefficients in salt mixtures [50]. MC simulations can be very useful for mixed salts because a salt mixture is easily created by just adding different types of ions in the simulation box and interactions between co-ions as well as between co-ions and counter-ions are included. However, care must be taken to check if the simulation has converged to the global energy minimum and is not trapped in a local minimum. In this work, convergence to the global energy minimum was checked very carefully by performing simulations of different MC configurations. A reliable result was accepted when two simulations of different configurations gave the same result. In Fig. 4 the calculated mean ionic activity coefficients of 0.01 mol·L−1 HCl at increasing NaCl concentration up to 4 mol·kg−1 are compared with experimental data [51]. The simulations were performed by using the concentration dependent dielectric constant of NaCl solutions [46]. The ionic radii of NaCl were the same as used for single salt solutions, i.e., the Cl− radius is 0.181 and Na+ radius is 0.202 nm, the radius of H+ was adjusted to get a good fit. The fitted proton radius was 0.27 nm, which is larger than in pure HCl, i.e., 0.23 nm. However, for pure HCl solutions the activity coefficients were fitted at fixed dielectric constant of water 78.39 at 25 °C. A nearly perfect fit was obtained in the whole concentration range, i.e., up to 1 mol·L−1, see Sect. 3.5.

Fig. 4
figure 4

Mean ionic activity coefficients of 0.01 mol·L−1 HCl at increasing concentrations of NaCl calculated by the MC simulations are compared with the experimentally determined data [51]. I denotes the ionic strength given in molality (mol·kg−1), blue line/spheres experimental data and red spheres MC simulations data. MC simulations were performed at 25 °C, with concentration dependent dielectric constants of sodium chloride solutions according to [46] and the ionic radii of H+, Na+ and Cl− were fixed at 0.27, 0.202 and 0.181 nm, respectively

Full size image

3.4 Temperature Dependent Activity Coefficients

Activity coefficients of electrolytes at high temperatures are very relevant for industrial processes. For example, in chlorate production the reaction temperature is between 70 and 80 °C. It is well established by experimental measurements that with increasing temperature the dielectric constant of water decreases [52]. A decrease in dielectric constant can have considerable impact on the ion–ion interactions. For example, ionic solutions, which are predominantly composed of dissociated ions at room temperature, can have considerable amounts of ion pairing at high temperatures due to the decreased shielding of electrostatic interactions. A simple way to include the temperature effects in the modelling of ionic activity coefficients is to use the experimentally measured temperature dependent dielectric constant of water in MC simulations. In Fig. 5 the activity coefficients calculated by MC simulations for NaCl at 50 and 70 °C are compared with the experimental data.

Fig. 5
figure 5

Mean ionic activity coefficients calculated by MC simulations using the temperature dependent dielectric constant of water are compared with experimental data. Exp denotes experimental data from [46] and the concentration (c) is given as square root of mol·L−1. The radii of Na+ and Cl− used in the MC simulations at 50 and 70 °C were fixed at 0.168 and 0.181 nm, respectively

Full size image

Note that in this case a constant dielectric constant was used irrespective of the salt concentration. The reason for that is the lack of experiments showing the influence of salt concentration on the dielectric constant at elevated temperatures.

A good agreement between the calculated and measured activity coefficients was obtained up to 0.8 mol·L−1 at both temperatures. At concentrations higher than 0.8 mol·L−1 the MC simulations overestimated the experimental data. One possible reason could be that at high temperatures the percentage of ions existing in solution as ion pairs also increases with increasing ionic strength. Consequently, the mean ionic activity coefficient will be lowered, which is essentially what the experimental data show. Since such effects are not included in the primitive model used in the MC simulations, it is reasonable that MC simulations overestimate the experimental data. Another explanation could be that the dielectric constant decreases with increasing concentration as observed at room temperature. Similar results are found for HCl at 50 °C, i.e., the MC simulations overestimated the activity coefficient at higher concentrations, Fig. S1.

3.5 Single Ion Activity Coefficients

For complex electrolyte systems, such as the chlorate process or in seawater, it would be beneficial to use activity coefficients for the individual ions to estimate the activities in the mixed electrolytes. In the recent literature there has been a lively debate regarding the measurement of single ion activities. Some scientists have claimed that it is possible to measure the single ion activity coefficients by ion selective electrodes while others have argued that it is impossible to determine single ion activity coefficients [18, 19, 24, 28, 30, 53]. Monte Carlo simulations provide a unique opportunity to shed some light on this controversial issue because single ion activities are calculated. The single ion activity is otherwise difficult to calculate from theory [54]. Moreover, simulations are less approximate than theories and usually theoretical results are first evaluated against computer simulations before comparing them with experimental data [55]. We shall start our exploration of this interesting issue by simulating the mean ionic activity and osmotic coefficients of HCl at 25 °C.

Mean ionic activity coefficients and osmotic coefficients obtained by taking an appropriate average of single ion activity coefficients are compared in Fig. 6 with the experimental data obtained by vapor pressure measurements [47] and by ion selective electrodes [22]. The MC data were converted from the molar to molal scale by using the solution density as described in detail in our previous study [9]. The molal concentration scale was chosen in order to keep the experimental data in its original concentration scale. A very good agreement between the calculated activity coefficients and measured values is apparent up to 1 mol·kg−1, whereas in the case of osmotic coefficient very good agreement up to 3.5 mol·kg−1 was obtained. Such a good agreement indicates that the ionic radii used for protons and chloride ions in the MC simulations are reasonable. A natural next question is whether single ion activities calculated by the MC simulations also match with the single ion activities measured by ion selective electrodes.

Fig. 6
figure 6

a Mean ionic activity and b osmotic coefficients calculated by MC simulations are compared with experimental data obtained by vapor pressure [47] as well as ion selective electrode measurements [22]. The concentration (m) is given as square root of mol·kg−1. The radii used in the MC simulations were fixed at 0.23 nm for the proton and 0.181 nm for Cl−, temperature at 25 °C and dielectric constant of water at 78.36

Full size image

In Fig. 7a and b the single ion activity coefficients of H+ and Cl− calculated by MC simulations and by Bates method [16] are compared with the experimental data obtained by ion selective electrodes.

Fig. 7
figure 7

Calculated single ion activity coefficients of a H+ and b Cl− are compared with experimental data from ion selective electrodes [22, 23]. The concentration (m) is given as square root of mol·kg−1. Data denoted by Bates are the calculated single ion activity coefficients from mean activity coefficients by using the Bates model [16]; c single ion activity coefficients for the Cl− ion calculated by using a Cl− radius equal to 0.0001 nm and H+ radius 0.23 nm are compared with experimental data as well as Cl− activity coefficients obtained by using the crystallographic radius for Cl−

Full size image

The striking features, which are apparent from these figures, are that the proton activity coefficient from MC is the lowest whereas the chloride ion activity coefficient is the highest compared with the experimental values [25]. In a recent study Quinones et al. [56] applied a modified Poisson-Boltzmann theory and obtained good agreement with our simulation data for single and mean activity coefficients of 1:1, 2:1 and 3:1 electrolytes. The upward bending behavior of experimentally determined activity or osmotic coefficients in sufficiently concentrated solutions is generally considered to be due to ion size effects. For example, in highly dissociated solutions such as HCl the co-ion repulsions get stronger as the number of ions in the same volume increases, which leads to the upward bending of the activity coefficients [2].

This matter was further investigated by performing MC simulations in which the proton radius was 0.23 nm whereas the chloride radius was fixed at 0.0001 nm, i.e., almost like an ion without size. The chloride activity coefficients calculated by this model are shown in Fig. 7c. Although these simulations have captured the trend of the experimental data of chloride activity coefficients, it also means that the chloride ion has no almost size, which seems unphysical. Note that the mean activity coefficients obtained by these MC simulations and ion selective electrodes are in agreement, Fig. 6a, because the measured chloride and proton activity coefficients are too low and too high, respectively, see Fig. 7a and b.

The intricate relation between proton and chloride has recently been highlighted for HCl by Reed et al. [57, 58]. The IR data were interpreted such that the Cl− ion resides inside the structure of a hydrated proton. This means that the effective size of a hydrated proton is much larger than what is normally considered, i.e., the proton size in aqueous solutions is close to the size of a hydrated lithium ion [5]. If this is the case, then the high measured values of proton activity coefficients and very low values of Cl− activities can be rationalized. Moreover, the presence of ion pairs in 2 mol·L−1 HCl has been reported [59] as well as ion-pairing in 0.1 mol·L−1 HCl at the water–atmosphere interphase [60]. The fact that the chloride ion in HCl has no effective size, as determined in the present study, might therefore not be totally “unphysical”.

Bates et al. [16, 61] provided a simple method by which single ion activity coefficients can be calculated from mean values if ionic hydration numbers are known. In this model the hydration number of chloride ion is considered to be zero, whereas for protons, sodium and potassium ions the hydration numbers are 8, 3.5 and 1.5, respectively. From Fig. 8a and b it is apparent that results based on the Bates model are close to the MC simulation data. In a recent study it was shown that the proton activity coefficients in HCl measured by using the ionic liquid based salt bridge are in good agreement with the prediction of the SIS model [23]. However, the SIS model could not predict the activity coefficient of protons in HBr measured by the same method [62]. The results for acids seem to be more complex than for the corresponding salts. In Fig. 8a and b the single ion activity coefficients of Na+ and Cl− ions calculated by MC simulations are compared with experimental data [21] and the Bates model.

Fig. 8
figure 8

Calculated single ion activity coefficients of a Na+ and b Cl− are compared with experimental data from ion selective electrodes [21]. The concentration (m) is given as square root of mol·kg−1. Data denoted by Bates are the calculated single ion activity coefficients from mean activity coefficients by using the Bates model. The radii used in the MC simulations were fixed at 0.168 nm for Na+ and 0.181 nm for Cl−, the temperature at 25 °C and dielectric constant of water at 78.36

Full size image

In this case the agreement between the experiments and MC data for the Na+ ion is much better than for H+ in HCl. It should be noted that no fitting procedure was used and the Na+ and Cl− radii used in the MC simulation are 0.168 nm and 0.181 nm, respectively. However, the experimentally determined chloride activity coefficients are still lower than the MC data. Interestingly, the trend in the chloride activity coefficients in KCl (Fig. S2b) is the opposite. In KCl the experimentally determined chloride activity coefficients are higher than the MC values but the K+ activities from the MC simulation and experiment are in reasonable agreement at least at low concentrations (Fig. S2a). These results indicate that there is an intricate coupling of the activities of ions in electrolytes. In a recent study [54] the nature of single ion activity coefficients was investigated by using MC simulations including ion–ion and ion–solvent interactions. The ion–solvent interactions were included by using the Born model and a concentration dependent dielectric constant of solutions was used. Reasonable agreement was found for 2:1 salts such as CaCl2 but very poor agreement for 1:1 salts such as NaCl. In our previous study [9] we have shown that the Born model is not good enough to include ion–solvent interactions.

3.6 Calculation of Water Activity

In order to explore further whether the best-fit ionic radii obtained by MC simulations yield reasonable thermodynamic results for properties other than activity and osmotic coefficients, the activity of water and pH were calculated as a function of electrolyte concentration. The activity of water has been measured in HCl as well as in NaCl solutions by the hygroscopic method [63]. Measured water activities are in good agreement with the water activities obtained from vapor pressure measurements [34, 47].

In this work the water activity (aw) was calculated by using the osmotic coefficients from the MC simulations according to Eq. 6,

$$a_{w} = \exp \left( { - \frac{vm\phi }{55.51}} \right)$$
(6)

where ν is the number of ions formed when one mole of salt is dissolved (2 for a 1:1 electrolyte), m is the salt molality, ϕ is the osmotic coefficient and 55.51 is the inverse of the molar mass of water in mol·kg−1.

The water activities calculated for HCl and NaCl are compared with the experimentally determined values in Fig. 9a and b, respectively.

Fig. 9
figure 9

Activities of water in a HCl, b NaCl solutions calculated by the MC simulations are compared with the experimentally determined data [63]. The concentration (m) is given as mol·kg−1. MC simulations were performed by using the ionic radii of H+, Na+, Cl−, 0.23, 0.168 and 0.181 nm, respectively; the temperature was fixed at 25 °C and the dielectric constant of water at 78.36

Full size image

An excellent agreement for HCl up to 3 mol·kg−1 concentration is obtained by using the ion radius of H+, i.e., 0.23 nm. This value of the fitted proton radius is in good agreement with the proton radii in different acids as was reported in our previous study [5]. In the case of NaCl very good agreement was found only up to 1.5 mol·kg−1 and at higher concentrations than that the calculated water activities are lower than the experimental ones. This best fitted range coincides with the range in which activity and osmotic coefficients of NaCl can be best modelled by using primitive models.

Interestingly, in a recent study [64], good agreement between the experimentally determined water activity in NaCl solutions and calculated by molecular dynamics simulations (MD) was also restricted to 1.5 mol·kg−1. However, a good agreement between the mean ionic activity coefficients calculated by the MD simulations and experimental data was only found for concentrations less than 0.8 mol·L−1. This is explained by the model used in the MD simulations where ion size parameters were obtained by fitting to the density of aqueous solutions [64, 65].

3.7 Calculation of pH in Concentrated Electrolyte Solutions

The pH is defined as − log{H+} and will also depend on water activity. Sinanayke [66] used Eq. 6 to calculate the pH in chloride rich solutions, where h is the hydration number and γ the activity coefficient of the proton.

$${\text{pH }} = {-}{ \log }_{ 10} \frac{{\left[ {{\text{H}}\left( {{\text{H}}_{2} {\text{O}}} \right)_{h}^{ + } } \right] \cdot \gamma_{{{\text{H}}\left( {{\text{H}}_{2} {\text{O}}} \right)_{h}^{ + } }} }}{{(a_{\text{w}} )^{h} }}$$
(7)

In Fig. 10 the pH values of HCl calculated by using Eq. 6 and the single ion activity coefficient of protons along with water activity values obtained by MC simulations are compared with the pH values measured by a hydrogen gas electrode [66, 67] as well as by a glass electrode [21]. The best fit to experimental pH data, measured both by glass electrode and hydrogen gas electrode, was obtained by including both hydration numbers (four for H+) and water activity in the MC simulations (Fig. 10). The results show that MC simulations can be used to estimate the pH in concentrated solutions, taking into account both hydration of the proton and the activity of water.

Fig. 10
figure 10

Comparison of the pH values in HCl solutions calculated by using the proton activity coefficients and water activity values obtained from MC simulations in Eq. 6 as given in [66] with the pH values determined by hydrogen gas and pH electrodes. MC simulations were performed at 25 °C and radii of the proton and chloride ions were fixed at 0.23 and 0.181 nm, respectively. Proton activity coefficients of Schneider are from Ref. [21] and of Senanayake are from Ref. [67]

Full size image

4 Conclusions

The applicability of the MC simulations method to calculate ionic activity, osmotic coefficient and water activity in concentrated salt solutions has been demonstrated, using only one adjustable parameter, i.e., the cation ion radius. Inclusion of a concentration dependent dielectric constant of water in the MC simulations was shown to extend its applicability range up to 3 mol·L−1 or beyond for some salts.

Temperature dependent activity coeffients could easily be calculated by using the temperature dependent dielectric constant of water while keeping only one adjustable parameter. The effective cation radii obtained by fitting the calculated data to the experimental data provide information about the state of ion hydration in solutions. It was shown that Li+ in LiClO3 solutions remains hydrated whereas Na+ in NaClO3 loses its hydration and strongly interacts with the chlorate ion.

The activity of water and pH in concentrated electrolytes were calculated from the MC simulations and excellent fits to experimental data were obtained. For pH both the water activity and the hydration of the proton were taken into account.

By comparing the single ion activity coefficients calculated by MC simulations, using the crystallographic radius for Cl− and 0.23 nm for the proton, with the experimentally determined values, it was shown that the measured individual activity coefficients for the proton are unrealistically high while for the chloride ion they are unrealistically low. A fit to the experimental data reveals that the chloride radius needs to be almost zero. This implies that the chloride ion is part of the hydration shell of the proton, a result that has been suggested in the literature.

The combination between experimental results and MC simulations is an efficient way to learn more about single activity coefficients and to shed some light on the controversy related to determine these values experimentally.