# Thermodynamic modeling of aqueous and mixed-solvent alkali chloride solutions using an ion-pairing equation of state

Abtin Raeispour Shirazi, Fufang Yang, Tri Dat Ngo, Nicolas Ferrando, Olivier Bernard, Jean-Pierre Simonin, Jean-Charles de Hemptinne

► To cite this version:

Abtin Raeispour Shirazi, Fufang Yang, Tri Dat Ngo, Nicolas Ferrando, Olivier Bernard, et al.. Thermo-dynamic modeling of aqueous and mixed-solvent alkali chloride solutions using an ion-pairing equation of state. Fluid Phase Equilibria, 2025, 588, pp.114231. 10.1016/j.fluid.2024.114231. hal-04941199

HAL Id: hal-04941199

https://ifp.hal.science/hal-04941199v1

Submitted on 11 Feb 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

# Thermodynamic modeling of aqueous and mixed-solvent alkali chloride solutions using an ion-pairing equation of state

Abtin Raeispour Shirazi<sup>a,b</sup>, Fufang Yang<sup>c</sup>, Tri Dat Ngo<sup>a</sup>, Nicolas Ferrando<sup>a</sup>, Olivier Bernard<sup>b</sup>, Jean-Pierre Simonin<sup>b</sup>, Jean-Charles de Hemptinne<sup>a,*</sup>

## ARTICLE INFO

**Keywords:**
Ion pairing
ePPC-SAFT
Binding MSA
Mixed-solvent
Wertheim theory
Bjerrum theory
RSP models

## ABSTRACT

In this work, the Binding-MSA (BiMSA) theory is implemented in an equation of state for the first time. The proposed model, BiMSA-electrolyte polar perturbed chain statistical associating fluid Theory (BiMSA-ePPC-SAFT), is applied to aqueous and mixed solvent electrolyte systems (water + methanol and water + ethanol) to investigate the impact of ion pairs. In a first step, a comparison is made between the Bjerrum and Wertheim theories for the calculation of the ion-ion association strength. The results obtained show that the Bjerrum theory is more successful in describing the association of ions especially in mixed solvent systems. In a second step, different types of relative static permittivity (RSP) models are implemented and compared. The obtained results reveal that using the Bjerrum theory with a volume-dependent RSP yields an ion-ion association strength that strongly changes with salinity. The models are further analyzed focusing on the relative importance of the various types of association (ion-solvent, solvent-solvent, solvent-cosolvent and cosolvent-cosolvent). It was observed that in an aqueous solution, ion-solvent and solvent-solvent association bonds are stronger and more important than those of ion pairing. However, for the mixed solvent systems, in high alcohol concentration, ion pairs exhibit the strongest bond.

## 1. Introduction

Thermodynamic modeling of electrolyte solutions stands as a cornerstone in various scientific and engineering disciplines, offering insights into the behavior of ionic solutions. It finds extensive applications across various fields, including underground storage in salt caverns [1], biorefining industry [2], pharmaceutical industry [3], wastewater treatment [4] and metal recovery using hydrometallurgy [5]. Central to the utility of thermodynamic modeling is its ability to accurately predict phase equilibria, which dictate the coexistence of different phases within a system. In electrolyte solutions, phase equilibria are affected by ion concentration, temperature, and pressure, demonstrating the complex interactions between electrostatic forces and ion-ion interactions. Of particular significance in many industrial and environmental contexts is the prediction of bubble pressure and vapour composition (vapour-liquid equilibria, VLE). These equilibria govern the distribution of species between the liquid and vapour phases, and their accurate prediction is paramount for the design of separation processes and the

understanding of environmental transport phenomena.

The mean ionic activity coefficient (MIAC) serves as a pivotal property in the characterization of electrolyte solutions, quantifying deviations from ideal behavior exhibited by ions in solution. This coefficient plays a vital role in elucidating ion interactions, solubility phenomena, and electrical conductivity [6], thus underpinning the design and optimization of electrochemical devices such as batteries and fuel cells. In addition to phase equilibria and activity coefficients, accurate prediction of density is crucial for optimizing mass transfer and transport processes. The density of electrolyte solutions is linked to factors such as ion size, charge, and solvent-solute interactions. Equation of state (EoS) models offer a systematic framework for the estimation of phase equilibria, VLE, MIAC, and density. EoS models encompass the thermodynamic behavior of complex systems, enabling the extrapolation of properties beyond experimental data and facilitating insights into system behavior under diverse conditions.

The thermodynamic modeling of electrolyte solutions has been extensively explored in academic literature, with notable reviews by Loehe and Donohue [7], Pinsky and Takano [8], Kontogeorgis et al.

<footnote><sup>a</sup> IFP Energies Nouvelles, 1 et 4 Avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France</footnote>

<footnote><sup>b</sup> CNRS, Laboratoire PHENIX, Campus P.M. Curie, Sorbonne Universités, F-75005 Paris, France</footnote>

<footnote><sup>c</sup> Laboratoire de Chimie de l'ENS Lyon, CNRS and Université de Lyon, 46 allée d'Italie, 69364 Lyon, France</footnote>

<footnote>* Corresponding author.
E-mail address: jean-charles.de-hemptinne@ifpen.fr (J.-C. de Hemptinne).</footnote>
Nomenclature

<table><tr><td>A</td><td>Helmholtz energy</td></tr><tr><td>a</td><td>activity</td></tr><tr><td>B</td><td>reducing coefficient in the Bjerrum association model</td></tr><tr><td>c</td><td>molarity</td></tr><tr><td>d<sub>v</sub></td><td>parameter in the Schreckenberg RSP model</td></tr><tr><td>d<sub>T</sub></td><td>parameter in the Schreckenberg RSP model</td></tr><tr><td>e</td><td>electronic charge</td></tr><tr><td>g<sup>hs</sup><sub>ij</sub></td><td>hard sphere i-j radial distribution function</td></tr><tr><td>K<sub>AB</sub></td><td>ion pair chemical equilibrium constant</td></tr><tr><td>k<sub>B</sub></td><td>Boltzmann constant</td></tr><tr><td>l<sub>B</sub></td><td>bjerrum length</td></tr><tr><td>m</td><td>molality</td></tr><tr><td>MIAC</td><td>mean ionic activity coefficient</td></tr><tr><td>N<sub>A</sub></td><td>Avogadro's number</td></tr><tr><td>OF</td><td>objective function</td></tr><tr><td>RSP</td><td>relative static permittivity</td></tr></table>

<table><tr><td>T</td><td>temperature</td></tr><tr><td>x</td><td>mole fraction</td></tr><tr><td>y</td><td>mole fraction in vapour phase</td></tr><tr><td>Z</td><td>valence</td></tr><tr><td>Z'</td><td>the effective charge of ions</td></tr><tr><td>X<sup>Ai</sup></td><td>fraction of non-bonded sites of A<sub>i</sub></td></tr></table>

Greek letters

<table><tr><td>γ</td><td>activity coefficient</td></tr><tr><td>Δ<sup>AiBj</sup></td><td>association strength</td></tr><tr><td>ε<sup>AB</sup></td><td>association energy parameter</td></tr><tr><td>ε<sub>0</sub></td><td>vacuum permittivity</td></tr><tr><td>ε</td><td>relative static permittivity</td></tr><tr><td>κ<sup>AB</sup></td><td>association volume parameter</td></tr><tr><td>ρ</td><td>density</td></tr><tr><td>σ<sup>hs</sup></td><td>hard sphere diameter</td></tr><tr><td>σ<sub>p</sub></td><td>sum of anion and cation diameters</td></tr></table>

[9–11], Ahmed et al. [12], Held [13], among others. Activity coefficient models and equations of state have been utilized to incorporate electrolyte terms, aiming to elucidate long-range ion-ion and ion-solvent interactions.

Some studies have employed non-primitive mean spherical approximation (MSA) models to explicitly consider these interactions [14–16], while others have relied on primitive MSA [17] or Debye-Hückel (DH) [18] theories within EoS frameworks, treating the solvent as a dielectric continuum. Additional terms, such as dispersion, association, Born, and multipolar terms, have been introduced to account for ion-dipole and dipole-dipole interactions [19–22]. Maribo-Mogensen et al. [23] conducted a comparative analysis between MSA and DH theories, noting their similarities within the EoS framework. Recently, another comparison between the two aforementioned theories was made by Simonin and Bernard [24] and the advantages of MSA theory in comparison to DH for the estimation of individual ion activity coefficients were discussed. The electrolyte term, either based on the Debye-Hückel theory or on the primitive MSA, is originally developed in a fixed volume framework (T, V, n), rather than fixed pressure (T, P, n). In what follows, we will focus on equation of state approaches that are also constructed as functions of volume (T, V, n).

Many electrolyte-EoSs have been proposed in the literature. Cubic plus association (CPA) + full version of the DH theory and Born term, was used by Kontogeorgis et al. [9,10,25] for the modeling of the solubility of salts in mixed solvents, VLE of water-methanol-LiCl, and liquid-liquid equilibria (LLE) of various mixed solvent systems (water-propanol-NaCl, water-propanol-n-octane-NaCl, water-1-butanol-KBr). In their works, in all cases ion-ion and ion-water interaction were described using the dispersion term and a composition dependent relative static permittivity (RSP) model. Das et al., [26,27] introduced the SAFT-VR + DE which is a nonprimitive electrolyte SAFT EoS. In their model a nonprimitive MSA version was employed to consider the ion-ion interactions. They used their model for water-salt systems up to very high concentrations and mixed solvent systems. SAFT-VRE, a combination of SAFT-VR, MSA and Born terms was used by Schreckenberg et al. [28] for the VLE and LLE modeling of various mixed solvent systems such as water-methanol-salt and water-1-butanol-salt. They introduced a new RSP model that depends on temperature, density and composition of the solvent. Selam et al. [29] introduced another extension for the SAFT-VR by adding DH and Born terms called eSAFT-VR Mie. In their model, the RSP is only solvent dependent. This model was also extended by Novak et al. [30] for calculation of gas solubilities in aqueous electrolyte solutions. Later on, they used their model for mixed solvent systems and introduced a

salt-concentration-dependent RSP [31,32]. Ahmed et al. [26] introduced an ePPC-SAFT model which is a combination of PC-SAFT, polar, MSA and Born terms. They used a composition dependent RSP and applied their model for mixed solvent systems water-alcohols-salt and CO₂-water-methanol-NaCl. Later on, their model was extended by Yang et al. [33,34] and applied to various mixed solvent systems. They introduced a new composition dependent RSP model, and the ion pairing phenomenon was discussed in their work. Held et al. [35] introduced ePC-SAFT which is a combination of PC-SAFT and DH term; their model was applied to mixed solvents systems with a constant RSP (only solvent concentration dependent). Later on, they introduced a new version of their previous model in which a Born term is added [36,37]. They considered the RSP as a salt-composition-dependent function which has a direct impact on the Born term. It was concluded that the Born term is efficient when the RSP is considered as a salt composition-dependent function.

It should be noted that in recent papers authors generally follow two different approaches: (i) one that considers dispersion and electrostatic interactions for ions and ignores any type of association [11,13,32,37]. No association is considered between ions or between ions and solvent molecules; (ii) one that ignores the dispersive interaction of ions and considers association (ion-ion and ion-solvent) and electrostatic interaction for the ions [33,34,38]. In recent papers by Bülow et al. [36], and Boroujeni et al. [39] the formation of ion pairs was discussed using the Bjerrum [40] model. The Bjerrum model was added into the EoS framework to consider electrostatic association of ions (ion pair formation). In their work, dispersive interactions between ions and ion-solvent and associative interaction between only ions are considered while ion-solvent association is ignored. Although considering the ion-solvent association is more realistic, one of the two interactions of ion-ion dispersion or ion-ion association could be ignored to reduce the number of adjustable parameters or complexity of the model. This was clearly discussed by Yang et al. [41]. They suggested that ignoring ion-ion association for the mixed solvents systems could be unreasonable or an oversimplification. The effect of ion-solvent association on the non-idealities in electrolyte systems was also pointed out by Wexler et al. [42].

The description of mixed solvent electrolytes has been identified in recent papers as very difficult, unless some empirical modification of the RSP is considered [11,32–34]. They may be taken as only temperature and solvent-composition dependent, they may be taken as volume dependent, or they may be taken with an explicit salt concentration dependence in order to improve the accuracy of their model with experimental data. Although introducing new RSP models increases the
<figure><img src="image_4.png" /><figcaption>Fig. 1. In the extended version of MSA theory (BiMSA), the effect of the formation of ion pairs and the long-range electrostatic interaction between ion pairs and free ions are added to the model.</figcaption></figure>

model accuracy, this should not be considered as the ultimate solution. We believe that there are still some missing pieces of this puzzle that require further investigation, and an extension of the electrolyte thermodynamic models other than through ‘effective’ RSP models (RSP models compensate the lack of different type of interactions in electrolyte models).

In parallel to equation of state development, it must be pointed out that the MSA theory itself has benefited from further developments, allowing it to consider more complex phenomena. This model was proposed initially by Percus and Yevick [43]. It was subsequently first solved by Waisman and Lebowitz in 1972 [44], and developed mainly by Blum [45]. It shares with the DH approximation the property of being a linearized theory, and of giving analytical expressions for the deviations from ideality in electrolyte solutions. These expressions are similar in form to the DH ones. However, the MSA was introduced as a refinement over the DH theory by accounting in a better way for the effect of the excluded volume of the ions on the correlation functions. At low concentration, it yields numerical results that are very close to those obtained from DH. At high concentration, the steric repulsion between the ions is expected to modify their distribution and the screening length. Later on, the MSA theory was extended by Bernard and Blum to account for the association of ions [46,47]. They solved the MSA equations for an arbitrary mixture of charged hard spheres with saturating bonds in the Wertheim [47] formalism when only pairs are allowed [48] and introduced the Binding version of the MSA (BiMSA). As a consequence of this discussion, one may point to several essential debates and concerns, namely:

* Which model (Wertheim or Bjerrum) best describes ion pairing in solution?

* What is the best RSP functionality?

* What type of interactions (ion-ion, ion-solvent, solvent-solvent, solvent-cosolvent or co-solvent-cosolvent) are more important for the modeling of aqueous and mixed solvent systems?

The focus of this work is to introduce the BiMSA theory to include ion pairs in a thermodynamic model for electrolyte solutions. To do so, mixed solvent systems comprising an alcohol as the co-solvent are modeled and discussed. As is well known, the presence of alcohols in the system drastically reduces the dielectric constant of the system and as a result, the concentration of ion pairs is expected to increase significantly [49,50].

In this work, the impact of ion pairs in mixed solvent solutions is investigated. To this end, a new framework is introduced for the ePPC-SAFT model, inspired by the Binding-MSA (BiMSA) theory introduced by Bernard et al. [46,51–53]. The proposed model is applied to chloride

salts-water-alcohol systems (these are the salts for which the most data are available). In the new framework, the formation of ion pairs using different association expressions such as Wertheim and Bjerrum is discussed and compared. The consequences of these choices on the RSP are discussed.

## 2. Framework

In the ePPC-SAFT model [33,34,38,41,54], the residual Helmholtz energy is expressed as a function of temperature, volume and mole numbers as follows:

$$A^{\text{res}} = A^{\text{hc}} + A^{\text{disp}} + A^{\text{assoc}} + A^{\text{polar}} + A^{\text{MSA}} + A^{\text{Born}} \quad (1)$$

where the superscripts hc, disp, assoc, and polar denote hard-chain, dispersion, association, and multipolar (as presented by Ref [21,55, 56]) effects, respectively, and Born describes the electrostatic contribution to the solvation [19]. The BiMSA superscript denotes the Binding MSA model [46,47]. This latter model was never combined with an equation of state which is therefore an originality of this work. It is explained and detailed in this section [51,52,57,58]. In the new framework, in order to introduce an ion pairing model, Aassoc and AMSA contributions within the ePPC-SAFT model are reconstructed.

### 2.1. Reconstruction of the MSA contribution ($A^{MSA}$)

In previous works [33,34], simplified MSA (SMSA) was used within the ePPC-SAFT framework and applied to various mixed-solvent systems. Despite a good accuracy of the model for mixed-solvent and aqueous systems, in the SMSA model electrostatic interactions are considered only between the free ions but not between the ion pairs and other ionic species. Thus, another version of the MSA theory was used, which explicitly includes an extension for the ion pairing systems by Bernard et al. [46–48]. This theory, called Binding MSA (BiMSA) was used to consider ion pair electrostatic interaction in the ePPC-SAFT framework. Therefore, as can be seen in Fig. 1, additional long-range interactions between the free ions and ion pairs are taken into account.

The BiMSA contribution was obtained from the Wertheim theory [59]. Here, similarly to SMSA as compared to the full MSA [60], we consider a ‘simplified’ version of BiMSA, in which the change in the Helmholtz free energy due to ion pairing is expressed as:

$$\frac{A^{\text{SBIMSA}}}{k_B T} = -\frac{e^2}{4\pi\epsilon_0 e k_B T} \left\{ \sum_{i=1}^{n_{\text{ion}}} \frac{n_i z_i^2 \Gamma_i}{1 + \Gamma_i \sigma_i} \right\} + \frac{V \Gamma_i^3}{3\pi} + \frac{A^{\text{MAL}}}{k_B T} \quad (2)$$

where AMAL term (Mass Action Law” contribution) represents the ion-ion association. While it belongs to the BiMSA theory, it is described further
<table><caption>**Table 1**
Temperature coefficients of dielectric constants in Eqs. (2), (3), (27) and (30) for the pure solvents.</caption><thead><tr><th>Solvent</th><th>d<sub>1</sub></th><th>d<sub>2</sub></th><th>d<sub>3</sub></th><th>d<sub>4</sub></th><th>d<sub>5</sub></th><th>T range(K)</th></tr></thead><tbody><tr><td>Water</td><td>-19.2905</td><td>29,814.5</td><td>-0.01968</td><td>0.013189</td><td>-0.03114</td><td>288.15-403.15</td></tr><tr><td>Ethanol</td><td>175.72</td><td>-3.0699</td><td>-0.3535</td><td>-0.20285</td><td>0.50644</td><td>288.15-328.15</td></tr></tbody></table>

<table><caption>**Table 2**
Types of interactions for the ion pairing approach in pure water solvent. '0' means no interaction, 'W' means the use of the Wertheim equation (Eq (18)) for computing the association strength, while 'B' means that the association strength is computed using the Bjerrum equation (Eq (26)). nc+ and na- are the number of solvating sites on each ion (interaction with solvents). In this work, nc+ = 7 and na- = 6, while ip=1, the ion pairing site number.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">water</th><th colspan="2">cation</th><th colspan="2">anion</th></tr><tr><th>A<br>++</th><th>B<br>-</th><th>C<br>nc<sup>+</sup></th><th>C'<br>ip<sup>+</sup></th><th>D<br>na<sup>-</sup></th><th>D'<br>ip<sup>-</sup></th></tr></thead><tbody><tr><td rowspan="2">water</td><td>A</td><td>++</td><td>0</td><td></td><td></td><td></td><td></td></tr><tr><td>B</td><td>-</td><td>W</td><td>0</td><td></td><td></td><td></td></tr><tr><td rowspan="2">cation</td><td>C</td><td>nc<sup>+</sup></td><td>0</td><td>W-solv</td><td>0</td><td></td><td></td></tr><tr><td>C'</td><td>ip<sup>+</sup></td><td>0</td><td>0</td><td>0</td><td></td><td></td></tr><tr><td rowspan="2">anion</td><td>D</td><td>na<sup>-</sup></td><td>W-solv</td><td>0</td><td>0</td><td>0</td><td></td></tr><tr><td>D'</td><td>ip<sup>-</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

when presenting the association (Section 2.3.2). The screening parameter is introduced using the following expression for a salt mixture:

$$4\Gamma^2 = \frac{e^2}{\epsilon_0 \epsilon k_B T} \sum_{i=1}^{n_{\text{ion}}} \frac{n_i}{V} \left( \frac{Z_i}{1 + \Gamma \sigma_i} \right)^2 + \frac{2e^2}{\epsilon_0 \epsilon k_B T V} \sum_{a_j=a_1}^{d_{\text{max}}} \times \sum_{c_k=\epsilon_1}^{\epsilon_{\text{max}}} \frac{n_{a_j c_k}}{\sigma_{p_{j,k}}} \left( \frac{\sigma_{a_j}}{1 + \Gamma \sigma_{a_j}} + \frac{\sigma_{c_k}}{1 + \Gamma \sigma_{c_k}} \right) Z'_{a_j} Z'_{c_k} \quad (3)$$

$$Z'_{a_j} = \frac{Z_{a_j}}{1 + \Gamma \sigma_{a_j}} \quad \text{and} \quad Z'_{c_k} = \frac{Z_{c_k}}{1 + \Gamma \sigma_{c_k}} \quad (4)$$

$$\sigma_{p_{j,k}} = \sigma_{a_j} + \sigma_{c_k} \quad (5)$$

Here, we limit ourselves to the study of a binary electrolyte. Then the screening parameter is determined using the following expression for a single salt:

$$4\Gamma^2 = \frac{e^2}{\epsilon_0 \epsilon k_B T} \sum_{i=1}^{n_{\text{ion}}} \frac{n_i}{V} \left( \frac{Z_i}{1 + \Gamma \sigma_i} \right)^2 + \frac{2e^2}{\epsilon_0 \epsilon k_B T V} \frac{n_{a_1 c_1}}{\sigma_{p_{1,1}}} \left( \frac{\sigma_{a_1}}{1 + \Gamma \sigma_{a_1}} + \frac{\sigma_{c_1}}{1 + \Gamma \sigma_{c_1}} \right) Z'_{a_1} Z'_{c_1} \quad (6)$$

where the subscript c denotes a cation, and the subscript a an anion. e is the electronic charge, ε is the relative permittivity of the medium, ε₀ the vacuum permittivity, Zᵢ is the valence, σᵢ is the diameter of the ion, nᵢ is the number of type i ions, Zₐⱼ and Zₐⱼ are the effective charge of ions (anions and cations respectively), σₚⱼₖ is the sum of anion and cation diameters and nₐⱼcₖ is the number of pairs formed between cations and anions introduced in Section 2.3.2. The screening parameter (Γ) is solved by iteration starting for the initial value Γ=κ/2 at infinite dilution (Γ→κ/2 at infinite dilution). In this case ions are treated as spheres of identical diameter σ. κ is the Debye screening length.

## 2.2. Relative static permittivity (dielectric constant)

As explained in the introduction, several investigations point to the key impact of the RSP functionality to the good representation of electrolyte systems: should it be salt-concentration dependent, or should it be volume (density) dependent; should it be solvent composition dependent? The need to consider the solvent composition dependence is settled in the literature. In this work two choices for RSP dependency are

considered. First, RSP is assumed to be volume and solvent composition dependent and second, it is solvent and temperature dependent but neither salt composition nor volume dependent. We have built our approach on the assumption discussed by several authors [11,32–34,38, 61] that a salt-concentration dependence needs not have the experimentally observed behaviour. The volume-dependence as used in other work [11,32–34,38,41,62,63] is considered because it would allow taking into account the continuity between liquid and vapour phases. In our work, two types of RSP models are used: Schreckenberg [28] which is volume-dependent and either Zuo & Fürst [64] (for water and ethanol) or Maryott & Smith [65] (for methanol) that implies that RSP is only temperature and solvent-composition dependent. In this work we never consider an explicit salt-concentration dependence.

### 2.2.1. Schreckenberg model

The Schreckenberg model [28] is a volume and solvent-composition dependent model which is expressed by the following equation:

$$\epsilon_r = 1 + \frac{n_{\text{solv}}}{V} d_v \left( \frac{d_T}{T} - 1 \right) \quad (7)$$

where n_solv is the number of moles of the solvent in volume V, d_v and d_T are parameters determined using the mixing rules.

$$d_v = \frac{\sum_{\text{solv}} n_i d_{v,i}}{\sum_{\text{solv}} n_i} \quad (8)$$

$$d_T = \frac{\sum_{\text{solv}} n_i d_{T,i}}{\sum_{\text{solv}} n_i} \quad (9)$$

where d_v and d_T values for water, methanol and ethanol are presented in supplementary information.

### 2.2.2. Zuo & Fürst and Maryott & Smith models

The Zuo & Fürst [64] model is only temperature-dependent and expressed by the following equation for the pure ethanol and water:

$$D_s = d_1 + \frac{d_2}{T} + d_3 T + d_4 T^2 + d_5 T^3 \quad (10)$$

where d₁-d₅ are constants that are fitted using experimental data [66, 67] (Table 1).

For pure methanol, the Maryott and Smith [65] model is used:

$$\log D_s = 1.514 - 0.00264(T - 298.15) \quad (11)$$

For the mixed solvents systems, the following mixing rules is used.

$$\epsilon_r = \frac{\sum_{\text{solv}} n_i D_i}{\sum_{\text{solv}} n_i} \quad (12)$$

## 2.3. Reconstruction of association contribution (A<sup>assoc</sup>)

For the calculation of A<sup>assoc</sup>, the Wertheim theory [59,68] is used:

$$\frac{A^{\text{assoc}}}{k_B T} = \sum_i n_i \sum_{A_i} \left[ \left( \ln X^{A_i} - \frac{X^{A_i}}{2} \right) + \frac{1}{2} M_i \right] \quad (13)$$

$$X^{A_i} = \left[ 1 + N_{Av} \rho \sum_{B_j} (X^{B_j} \Delta^{A_i B_j}) \right]^{-1} \quad (14)$$
<table><caption>**Table 3**
Summary of model combinations.</caption><thead><tr><th colspan="2">Contributions:</th></tr></thead><tbody><tr><td>Hard sphere (hc):</td><td>PC-SAFT</td></tr><tr><td>Dispersion (disp):</td><td>PC-SAFT</td></tr><tr><td>Electrostatic solvation (Born):</td><td>Born</td></tr><tr><td>Electrostatic ion-ion interaction (MSA):</td><td>sBiMSA (Eq.2)</td></tr><tr><td>Association (Assoc) (ion-solvent and solvent-solvent):</td><td>Wertheim (Eq.18)</td></tr><tr><td>Mass action law (MAL) (ion-ion association):</td><td>Werheim (Eq.18)</td></tr><tr><td rowspan="2">RSP:</td><td>Schreckenberg</td></tr><tr><td>Bjerrum (Eq.28) (g<sup>el</sup>=1)</td></tr><tr><td rowspan="2">Model name:</td><td>Schreckenberg</td></tr><tr><td>Bjerrum (Eq.28) (g<sup>el</sup>=1)</td></tr><tr><td></td><td>Zuo & Fürst & Maryott & Smith</td></tr><tr><td></td><td>BT</td></tr><tr><td></td><td>Bjerrum * g<sup>el</sup> (eq.29)</td></tr><tr><td></td><td>Schreckenberg</td></tr><tr><td></td><td>BgS</td></tr><tr><td></td><td>BgT</td></tr></tbody></table>

<table><caption>**Table 4**
Number of datasets, number of data points, temperature range and maximum ion concentration (mole fraction) for each property of each aqueous electrolyte solution used in the parameter regression in this work. The references can be found in previous works [38].</caption><thead><tr><th>Water +</th><th>Property</th><th>No. datasets</th><th>No. data points</th><th>T (K)</th><th>Max. x<sub>ion</sub></th></tr></thead><tbody><tr><td rowspan="3">NaCl</td><td>MIAC</td><td>20</td><td>344</td><td>273.15–333.15</td><td>0.089</td></tr><tr><td>VLE p<sup>*</sup></td><td>3</td><td>166</td><td>294.34–368.07</td><td>0.091</td></tr><tr><td>Density</td><td>129</td><td>2578</td><td>253.00–473.15</td><td>0.098</td></tr><tr><td rowspan="3">KCl</td><td>MIAC</td><td>4</td><td>43</td><td>273.15–298.15</td><td>0.056</td></tr><tr><td>VLE p<sup>*</sup></td><td>8</td><td>192</td><td>298.15–374.15</td><td>0.095</td></tr><tr><td>Density</td><td>104</td><td>1615</td><td>273.15–473.15</td><td>0.11</td></tr><tr><td rowspan="3">RbCl</td><td>MIAC</td><td>1</td><td>21</td><td>298.15</td><td>0.076</td></tr><tr><td>VLE p<sup>*</sup></td><td>3</td><td>65</td><td>303.15–368.15</td><td>0.10</td></tr><tr><td>Density</td><td>14</td><td>168</td><td>288.15–328.15</td><td>0.12</td></tr><tr><td rowspan="3">CsCl</td><td>MIAC</td><td>4</td><td>51</td><td>298.15</td><td>0.049</td></tr><tr><td>VLE p<sup>*</sup></td><td>2</td><td>75</td><td>303.15–368.15</td><td>0.12</td></tr><tr><td>Density</td><td>30</td><td>351</td><td>273.85–473.15</td><td>0.19</td></tr></tbody></table>

<table><caption>**Table 5**
Summary of the mixed-solvent electrolyte properties (MIAC, VLE<sub>p,y</sub>*, and density data) included in the parameterization. Number of data points, temperature range and maximum ion and alcohol concentration (mole fraction) for each property of each mixed solvent electrolyte solution used in the parameter regression in this work. The references can be found in previous work [33,34].</caption><thead><tr><th>System</th><th>Property</th><th>No. data points</th><th>T (K)</th><th>Max. x<sub>ion</sub></th><th>Max. x<sub>alcohol</sub></th></tr></thead><tbody><tr><td rowspan="2">Water + Methanol + NaCl</td><td>MIAC</td><td>327</td><td>298.15–318.15</td><td>0.068</td><td>0.834</td></tr><tr><td>VLE p<sub>o</sub>, y<sup>*</sup></td><td>49</td><td>298.34–347.15</td><td>0.092</td><td>0.909</td></tr><tr><td rowspan="3">KCl</td><td>Density</td><td>173</td><td>288.15–328.15</td><td>0.114</td><td>0.905</td></tr><tr><td>MIAC</td><td>65</td><td>298.15</td><td>0.036</td><td>0.457</td></tr><tr><td>VLE p<sub>o</sub>, y<sup>*</sup></td><td>18</td><td>298.15–372.15</td><td>0.095</td><td>0.897</td></tr><tr><td rowspan="2">RbCl</td><td>Density</td><td>20</td><td>273.15–473.15</td><td>0.127</td><td>0.981</td></tr><tr><td>MIAC</td><td>76</td><td>298.15</td><td>0.068</td><td>0.691</td></tr><tr><td rowspan="2">Water + Ethanol + NaCl</td><td>MIAC</td><td>250</td><td>288.15–318.15</td><td>0.052</td><td>0.776</td></tr><tr><td>VLE p<sub>o</sub>, y<sup>*</sup></td><td>51</td><td>298.15–332.15</td><td>0.071</td><td>0.726</td></tr><tr><td rowspan="3">KCl</td><td>Density</td><td>70</td><td>298.15–313.15</td><td>0.082</td><td>0.779</td></tr><tr><td>MIAC</td><td>25</td><td>298.15</td><td>0.028</td><td>0.201</td></tr><tr><td>VLE p<sub>o</sub>, y<sup>*</sup></td><td>24</td><td>298.15–352.15</td><td>0.026</td><td>0.481</td></tr><tr><td rowspan="3">CsCl</td><td>Density</td><td>69</td><td>298.15–313.15</td><td>0.072</td><td>0.779</td></tr><tr><td>MIAC</td><td>56</td><td>298.15</td><td>0.080</td><td>0.475</td></tr><tr><td>Density</td><td>50</td><td>298.15–308.15</td><td>~0.870</td><td>0.089</td></tr></tbody></table>

<table><caption>**Table 6**
Ion-specific parameterization procedure of WS model using IP-Wertheim, and the rest (BS, BT, BgS and BgT) using IP-Bjerrum for the aqueous chloride solutions (NaCl, KCl, RbCl and CsCl).</caption><thead><tr><th rowspan="2">Models:</th><th rowspan="2">Association model</th><th colspan="2">Parameters</th></tr><tr><th>Pseudo-unary</th><th>Binary</th></tr></thead><tbody><tr><td>WS</td><td>IP-Wertheim</td><td>ε<sup>AB</sup><sub>Na<sup>+</sup></sub>, ε<sup>AB</sup><sub>K<sup>+</sup></sub>, ε<sup>AB</sup><sub>Rb<sup>+</sup></sub>, ε<sup>AB</sup><sub>Cs<sup>+</sup></sub>, ε<sup>AB</sup><sub>Cl<sup>+</sup></sub>, r<sup>MSA</sup><sub>Na<sup>+</sup></sub>, r<sup>MSA</sup><sub>K<sup>+</sup></sub>, r<sup>MSA</sup><sub>Rb<sup>+</sup></sub>, r<sup>MSA</sup><sub>Cs<sup>+</sup></sub>, r<sup>MSA</sup><sub>Cl<sup>+</sup></sub>, r<sup>HS</sup></td><td>w<sub>Na<sup>+</sup></sub>–Cl<sup>+</sup>, w<sub>K<sup>+</sup></sub>–Cl<sup>+</sup>, w<sub>Rb<sup>+</sup></sub>–Cl<sup>+</sup>, w<sub>Cs<sup>+</sup></sub>–Cl<sup>+</sup></td></tr><tr><td>BS, BT, BgS and BgT</td><td>IP-Bjerrum</td><td>ε<sup>AB</sup><sub>Na<sup>+</sup></sub>, ε<sup>AB</sup><sub>K<sup>+</sup></sub>, ε<sup>AB</sup><sub>Rb<sup>+</sup></sub>, ε<sup>AB</sup><sub>Cs<sup>+</sup></sub>, ε<sup>AB</sup><sub>Cl<sup>+</sup></sub>, r<sup>MSA</sup><sub>Na<sup>+</sup></sub>, r<sup>MSA</sup><sub>K<sup>+</sup></sub>, r<sup>MSA</sup><sub>Rb<sup>+</sup></sub>, r<sup>MSA</sup><sub>Cs<sup>+</sup></sub>, r<sup>MSA</sup><sub>Cl<sup>+</sup></sub>, r<sup>HS</sup></td><td>b<sub>Na<sup>+</sup></sub>–Cl<sup>+</sup>, b<sub>K<sup>+</sup></sub>–Cl<sup>+</sup>, b<sub>Rb<sup>+</sup></sub>–Cl<sup>+</sup>, b<sub>Cs<sup>+</sup></sub>–Cl<sup>+</sup></td></tr></tbody></table>

<table><caption>**Table 7**
Ion-specific pseudo-unary parameters of four models introduced in Table 3 for chloride salts in aqueous solution (the regressed parameters are marked as bold and italic).</caption><thead><tr><th>N<sub>site</sub></th><th>σ<sub>Pauling</sub> (Å)<sup>a</sup></th><th>r<sup>HS</sup></th><th>ε<sup>AB</sup>/k (K)</th><th>r<sup>MSA</sup></th><th>σ<sub>Born</sub> (Å)<sup>c</sup></th></tr></thead><tbody><tr><td colspan="6">WS model</td></tr><tr><td>Na<sup>+</sup></td><td>7</td><td>1.9</td><td><strong>0.955</strong></td><td><strong>1117.7</strong></td><td><strong>1.89</strong></td><td>3.235</td></tr><tr><td>K<sup>+</sup></td><td>7</td><td>2.66</td><td><strong>0.955</strong></td><td><strong>592.73</strong></td><td><strong>1.35</strong></td><td>3.897</td></tr><tr><td>Rb<sup>+</sup></td><td>7</td><td>2.96</td><td><strong>0.955</strong></td><td><strong>161.6</strong></td><td><strong>1</strong></td><td>4.169</td></tr><tr><td>Cs<sup>+</sup></td><td>7</td><td>3.38</td><td><strong>0.955</strong></td><td><strong>-261.7</strong><sup>d</sup></td><td><strong>2.5</strong></td><td>4.482</td></tr><tr><td>Cl<sup>-</sup></td><td>6</td><td>3.62</td><td><strong>0.955</strong></td><td><strong>1987</strong></td><td><strong>1.57</strong></td><td>4.512</td></tr><tr><td colspan="6">BS model</td></tr><tr><td>Na<sup>+</sup></td><td>7</td><td>1.9</td><td><strong>0.891</strong></td><td><strong>2546</strong></td><td><strong>2.74</strong></td><td>3.235</td></tr><tr><td>K<sup>+</sup></td><td>7</td><td>2.66</td><td><strong>0.891</strong></td><td><strong>1805</strong></td><td><strong>1.807</strong></td><td>3.897</td></tr><tr><td>Rb<sup>+</sup></td><td>7</td><td>2.96</td><td><strong>0.891</strong></td><td><strong>2501</strong></td><td><strong>1.28</strong></td><td>4.169</td></tr><tr><td>Cs<sup>+</sup></td><td>7</td><td>3.38</td><td><strong>0.891</strong></td><td><strong>-191</strong><sup>d</sup></td><td><strong>3.15</strong></td><td>4.482</td></tr><tr><td>Cl<sup>-</sup></td><td>6</td><td>3.62</td><td><strong>0.891</strong></td><td><strong>3237</strong></td><td><strong>1.014</strong></td><td>4.512</td></tr><tr><td colspan="6">BT model</td></tr><tr><td>Na<sup>+</sup></td><td>7</td><td>1.9</td><td><strong>0.937</strong></td><td><strong>1373</strong></td><td><strong>2.86</strong></td><td>3.235</td></tr><tr><td>K<sup>+</sup></td><td>7</td><td>2.66</td><td><strong>0.937</strong></td><td><strong>949</strong></td><td><strong>1.46</strong></td><td>3.897</td></tr><tr><td>Rb<sup>+</sup></td><td>7</td><td>2.96</td><td><strong>0.937</strong></td><td><strong>493</strong></td><td><strong>1.89</strong></td><td>4.169</td></tr><tr><td>Cs<sup>+</sup></td><td>7</td><td>3.38</td><td><strong>0.937</strong></td><td><strong>-155</strong><sup>d</sup></td><td><strong>2.33</strong></td><td>4.482</td></tr><tr><td>Cl<sup>-</sup></td><td>6</td><td>3.62</td><td><strong>0.937</strong></td><td><strong>2766</strong></td><td><strong>1.48</strong></td><td>4.512</td></tr><tr><td colspan="6">BgT model</td></tr><tr><td>Na<sup>+</sup></td><td>7</td><td>1.9</td><td><strong>0.952</strong></td><td><strong>1884</strong></td><td><strong>3.24</strong></td><td>3.235</td></tr><tr><td>K<sup>+</sup></td><td>7</td><td>2.66</td><td><strong>0.952</strong></td><td><strong>948</strong></td><td><strong>1.3</strong></td><td>3.897</td></tr><tr><td>Rb<sup>+</sup></td><td>7</td><td>2.96</td><td><strong>0.952</strong></td><td><strong>1620</strong></td><td><strong>2.15</strong></td><td>4.169</td></tr><tr><td>Cs<sup>+</sup></td><td>7</td><td>3.38</td><td><strong>0.952</strong></td><td><strong>1482</strong></td><td><strong>2.27</strong></td><td>4.482</td></tr><tr><td>Cl<sup>-</sup></td><td>6</td><td>3.62</td><td><strong>0.952</strong></td><td><strong>1928</strong></td><td><strong>1.3</strong></td><td>4.512</td></tr></tbody></table>

where M<sub>i</sub> stand for the number of association sites on molecule i, X<sup>Ai</sup> is the fraction of non-bonded sites of A<sub>i</sub> and Δ<sup>Ai,Bj</sup> is called the association strength and is in fact a density-based association constant (see below).

<footnote>* VLE *p* is the vapour pressure.</footnote>

<footnote>*VLE *p*, *y* is the vapour pressure and composition.</footnote>

<footnote>a The Pauling diameter values are obtained from Ref [74].</footnote>

<footnote>c The Born diameters are calculated using the Gibbs energy of solvation [76] using the procedure introduced in previous work [61].</footnote>

<footnote>d The pseudo-pure association energy parameter, ε<sup>AB</sup><sub>Cs<sup>+</sup></sub>, can be negative. However, ε<sup>AB</sup><sub>Cs<sup>+</sup>–w</sub> calculated using Eq. (33) is always positive.</footnote>
<table><caption>**Table 8**
Cation-anion association parameters for the four models introduced in Table 3 for chloride salts in aqueous solutions.</caption><thead><tr><th></th><th>WS</th><th>BS</th><th>BT</th><th>BgT</th></tr><tr><th></th><th>ε<sup>AB</sup><sub>c-a</sub>/k (K)</th><th>b<sub>c-cl</sub></th><th>b<sub>c-cl</sub></th><th>b<sub>c-cl</sub></th></tr></thead><tbody><tr><td>Na<sup>+</sup></td><td>1106.2</td><td>0.73</td><td>0.74</td><td>0.75</td></tr><tr><td>K<sup>+</sup></td><td>1257.5</td><td>0.82</td><td>0.84</td><td>0.86</td></tr><tr><td>Rb<sup>+</sup></td><td>1347.4</td><td>0.86</td><td>0.94</td><td>0.96</td></tr><tr><td>Cs<sup>+</sup></td><td>1635.6</td><td>1.15</td><td>1.12</td><td>0.95</td></tr></tbody></table>

### 2.3.1. Association scheme

In this work, we aim at describing both solvation and ion pairing using the association term (Eq. (1)). This is done by adapting the association scheme as follows:

As can be seen in Table 2, water is described using the 4C association scheme: two positive and two negative sites, which means that water can associate with any solvent molecule (e.g., water or alcohols) or with ions (negative sites of water with cations and positive sites of water with anions). In Table 2, where only water is considered as the solvent, these sites are labelled as A and B: two positive (electron acceptor) sites A and two negative (electron-donor) sites B. For anions, only negative sites and for cations only positive sites are considered which means that ions can associate with both solvent molecules and their counterions (ion-solvent and ion-ion association). Cations can thus form bonds with anions in an indefinite way, leading to ion clustering. This is something that can sometimes be observed but remains limited in practice [34]. In order to avoid ion-clustering, two different types of sites on the ions are distinguished: one type of site associates with unlike ions (they carry a prime in Table 2. In our work, only one type of such site exists on every ion: ip⁺ = ip⁻ = 1), and the other one associates with solvent molecules. In case of non-stoichiometric salts (e.g. 2–1) then either ip⁺ or ip⁻ could take another integer value depending on the number of ions that it can interact with. The number of ionic sites that interact with solvent

molecules are here indicated by nc⁺ and na⁻, respectively. In this work, for cations nc⁺ = 7 and for anions na⁻ = 6. These numbers are taken from the previous works [38,41] in our group and are in accordance with molecular simulation [69] or X-ray diffraction studies. Table 2 shows how these different types of sites can interact (it is a symmetric table so only half of the table is shown). When a zero is shown, no interaction occurs (or ΔA,Bi=0). The label ‘W’ means the usual hydrogen-bond for water molecules that is described with the regular Wertheim expression (Eq. (18)). Ion-solvent interactions are also expressed using the Wertheim Eq. (18). They are labelled as ‘W-solv’ as they describe the structure-breaking effect of solvation. Finally, the ion-ion interactions (‘B-IP’) are expressed using a different association strength (Bjerrum, Eq. (26)), as explained below.

As a conclusion, one can consider that the association term is constructed out of three contributions:

$$A^{\text{assoc}} = A^{\text{ion-ion}} + A^{\text{ion-solvent}} + A^{\text{solvent-solvent}} \quad (15)$$

The ion-solvent and the solvent-solvent contributions cannot be split in Eq. (15), because the same site on a solvent molecule can associate with either a solvent site or an ion site. Yet, the ion-ion sites (C' and D') only interact with other ion-ion sites, hence their contribution can be separated, which is called ion-ion association A<sup>ion-ion</sup>.

### 2.3.2. Association of ion pairs ($A^{MAL} = A^{ion-ion}$)

$A^{ion-ion}$ in equation 15 is exactly equal to $A^{MAL}$ contribution which comes from BiMSA model. In the BiMSA publication [51], the expression is:

$$\frac{A^{\text{MAL}}}{k_B T} = \sum_{i=1}^{n_{\text{ion}}} n_i \ln X^{A_i} + n_{a_1 c_1} \quad (16)$$

where $n_{a_1 c_1}$ is the number of pairs obtained from the following expression:

<table><caption>**Table 9**
Alcohol-ion energy parameters for the four models of Table 3 for the mixed-solvent electrolyte solutions.</caption><thead><tr><th rowspan="2">BIP</th><th colspan="4">Methanol</th><th colspan="4">Ethanol</th></tr><tr><th>WS</th><th>BS</th><th>BT</th><th>BgT</th><th>WS</th><th>BS</th><th>BT</th><th>BgT</th></tr></thead><tbody><tr><td>W<sub>Na⁺</sub>-alcohol</td><td>-0.122</td><td>0.277</td><td>0.013</td><td>-0.028</td><td>0.151</td><td>0.111</td><td>0.116</td><td>0.125</td></tr><tr><td>W<sub>K⁺</sub>-alcohol</td><td>-0.028</td><td>0.998</td><td>0.156</td><td>0.030</td><td>0.136</td><td>0.132</td><td>0.537</td><td>0.135</td></tr><tr><td>W<sub>Cl⁻</sub>-alcohol</td><td>0.161</td><td>-0.092</td><td>0.120</td><td>0.076</td><td>0.139</td><td>0.076</td><td>0.118</td><td>0.059</td></tr><tr><td>W<sub>Rb⁺</sub>-alcohol</td><td>-0.536</td><td>0.235</td><td>-0.071</td><td>0.110</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>W<sub>Cs⁺</sub>-alcohol</td><td>—</td><td>—</td><td>—</td><td>—</td><td>-0.071</td><td>0.710</td><td>-0.092</td><td>0.091</td></tr></tbody></table>

**Table 10**

Average absolute percentage deviations (AADs*) and maximum absolute percentage deviations (MADs*) for MIAC, VLE vapour pressure (p), and density of WS, BS, BT and BgT models for chloride salts in aqueous solutions.

<table><thead><tr><th rowspan="2">Model name</th><th colspan="4">MIAC</th><th colspan="4">VLE p</th><th colspan="4">Density</th></tr><tr><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th></tr></thead><tbody><tr><td>Water +</td><td>AAD(%)</td><td>MAD(%)</td><td>AAD(%)</td><td>MAD(%)</td><td>AAD(%)</td><td>MAD(%)</td><td>AAD(%)</td><td>MAD(%)</td><td>AAD(%)</td><td>MAD(%)</td><td>AAD(%)</td><td>MAD(%)</td></tr><tr><td>NaCl</td><td>0.35</td><td>2.84</td><td>1.21</td><td>8.56</td><td>1.2</td><td>2.27</td><td>0.93</td><td>2.25</td><td>0.54</td><td>3.83</td><td>0.52</td><td>2.56</td></tr><tr><td>KCl</td><td>0.43</td><td>1.72</td><td>2.48</td><td>7.81</td><td>1.63</td><td>4.35</td><td>1.39</td><td>3.22</td><td>0.22</td><td>1.65</td><td>0.94</td><td>4.25</td></tr><tr><td>RbCl</td><td>0.18</td><td>0.62</td><td>1.17</td><td>2.54</td><td>1.07</td><td>4.52</td><td>1.24</td><td>3.92</td><td>0.15</td><td>0.54</td><td>1.48</td><td>6.94</td></tr><tr><td>CsCl</td><td>0.6</td><td>1.35</td><td>0.49</td><td>1.35</td><td>1.09</td><td>6.61</td><td>0.75</td><td>2.88</td><td>0.29</td><td>5.05</td><td>1.05</td><td>7.69</td></tr><tr><td>Model name</td><td>BT</td><td>BgT</td><td>BT</td><td>BgT</td><td>BT</td><td>BgT</td><td>BT</td><td>BgT</td><td>BT</td><td>BgT</td><td>BT</td><td>BgT</td></tr><tr><td>Water+</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td></tr><tr><td>NaCl</td><td>1.1</td><td>7.84</td><td>1.14</td><td>6.49</td><td>0.88</td><td>2.02</td><td>0.93</td><td>2.15</td><td>0.48</td><td>3.86</td><td>0.56</td><td>4.56</td></tr><tr><td>KCl</td><td>0.42</td><td>1.48</td><td>1.05</td><td>2.1</td><td>1.41</td><td>4.28</td><td>1.36</td><td>3.24</td><td>0.3</td><td>1.22</td><td>0.27</td><td>2.37</td></tr><tr><td>RbCl</td><td>1.16</td><td>3.94</td><td>1.27</td><td>2.69</td><td>0.85</td><td>2.02</td><td>0.96</td><td>2.03</td><td>0.34</td><td>1.58</td><td>0.47</td><td>2.29</td></tr><tr><td>CsCl</td><td>0.6</td><td>1.22</td><td>1.23</td><td>6.63</td><td>0.86</td><td>3.41</td><td>1.11</td><td>3.23</td><td>0.45</td><td>12.5</td><td>0.83</td><td>3.97</td></tr></tbody></table>

$$ * \text{AAD}\% = \frac{\sum_{i=1}^{n_{dp}} \left| \frac{\Psi_{i, \text{cal}} - \Psi_{i, \text{exp}}}{\Psi_{i, \text{exp}}} \right|}{n_{dp}^{\Psi}} * 100, \quad \text{MAD}\% = \max \left| \frac{\Psi_{i, \text{cal}} - \Psi_{i, \text{exp}}}{\Psi_{i, \text{exp}}} \right| * 100, \quad \psi \text{ is MIAC, density or vapour pressure (VLE p).} $$
<caption>**Table 11**</caption>

Average absolute percentage deviations (AADs*) and maximum absolute percentage deviations (MADs*) for MIAC, VLE, and density of WS, BS, BT and BgT models for chloride salts in water+Ethanol and water+methanol systems. The MIAC, VLE *p*, and density deviations are in percentage. The VLE *y* absolute deviations are absolute deviations multiplied by 100.

<table><thead><tr><th rowspan="2">Model name</th><th colspan="4">MIAC</th><th colspan="4">VLE p</th><th colspan="4">VLE y</th><th colspan="4">Density</th></tr><tr><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th><th>WS</th><th>BS</th></tr></thead><tbody><tr><td>Water+<br>Ethanol +<br>NaCl</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td></tr><tr><td>KCl</td><td>14.78</td><td>49.12</td><td>6.24</td><td>21.31</td><td>2.07</td><td>10.96</td><td>2.63</td><td>14.45</td><td>2.46</td><td>5.06</td><td>2.67</td><td>6.49</td><td>0.83</td><td>2.33</td><td>0.74</td></tr><tr><td>CsCl</td><td>10.63</td><td>30.48</td><td>4.22</td><td>6.53</td><td>2.3</td><td>4.46</td><td>1.54</td><td>4.38</td><td>2.56</td><td>6.52</td><td>3.65</td><td>6.31</td><td>0.52</td><td>6.14</td><td>0.96</td></tr><tr><td></td><td>20.89</td><td>72.9</td><td>6.98</td><td>28.41</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.55</td><td>3.09</td><td>0.4</td></tr><tr><td>BT</td><td colspan="4">BdT</td><td colspan="4">BdT</td><td colspan="4">BdT</td><td colspan="4">BdT</td></tr><tr><td>AAD (%)</td><td>33.07</td><td>17.48</td><td>56.49</td><td>2.5</td><td>15.08</td><td>3.02</td><td>18.48</td><td>2.68</td><td>8.51</td><td>5.06</td><td>20.76</td><td>0.76</td><td>2.27</td><td>0.85</td><td>2.29</td></tr><tr><td>NaCl</td><td>12.19</td><td>20.192</td><td>11.44</td><td>28.7</td><td>1.23</td><td>4.8</td><td>2.77</td><td>5.7</td><td>3.34</td><td>5.5</td><td>4.19</td><td>10.16</td><td>0.56</td><td>6.57</td><td>0.61</td><td>5.83</td></tr><tr><td>KCl</td><td>8.89</td><td>20.192</td><td>11.44</td><td>28.7</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.58</td><td>3.18</td><td>0.29</td><td>2.69</td></tr><tr><td>CsCl</td><td>4.50</td><td>15.37</td><td>26.32</td><td>92.45</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Model name</td><td>WS</td><td>BS</td><td>WS</td><td>BS</td><td>WS</td><td>BS</td><td>WS</td><td>BS</td><td>WS</td><td>BS</td><td>WS</td><td>BS</td><td>WS</td><td>BS</td><td>WS</td><td>BS</td></tr><tr><td>Water+<br>Methanol+<br>NaCl</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td><td>AAD (%)</td><td>MAD (%)</td></tr><tr><td>KCl</td><td>6.58</td><td>31.88</td><td>3.94</td><td>18.54</td><td>2.68</td><td>13.58</td><td>2.37</td><td>10.41</td><td>1.07</td><td>5.06</td><td>1.23</td><td>4.78</td><td>2.68</td><td>13.58</td><td>1.9</td><td>8.87</td></tr><tr><td>RbCl</td><td>6.48</td><td>34.89</td><td>2.08</td><td>13.64</td><td>2.97</td><td>9.85</td><td>2.68</td><td>6.8</td><td>0.97</td><td>5.11</td><td>0.84</td><td>4.22</td><td>2.97</td><td>9.85</td><td>6.05</td><td>13.76</td></tr><tr><td></td><td>6.12</td><td>61.3</td><td>3.62</td><td>21.78</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>BT</td><td colspan="4">BdT</td><td colspan="4">BdT</td><td colspan="4">BdT</td><td colspan="4">BdT</td></tr><tr><td>AAD (%)</td><td>11.37</td><td>5.4</td><td>19.3</td><td>2.31</td><td>12.2</td><td>2.11</td><td>7.90</td><td>1.11</td><td>4.63</td><td>0.99</td><td>4.07</td><td>1.3</td><td>4.7</td><td>1.26</td><td>4</td><td></td></tr><tr><td>NaCl</td><td>3.54</td><td>9.56</td><td>5.03</td><td>19.75</td><td>3.10</td><td>8.73</td><td>2.54</td><td>8.47</td><td>1.12</td><td>5.12</td><td>2.12</td><td>9.03</td><td>3.8</td><td>9.04</td><td>2.97</td><td>6.51</td></tr><tr><td>KCl</td><td>3.23</td><td>9.56</td><td>5.03</td><td>19.75</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>RbCl</td><td>2.45</td><td>8.49</td><td>3.55</td><td>18.89</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

* For the vapour composition (VLE *y*): AAD% = $\frac{\sum_{i=1}^{n_y^p} |y_i \text{ cal} - y_i \text{ exp}|}{n_{dp}^y} * 100$, MAD% = $\max|y \text{ cal} - y \text{ exp}| * 100$ and for the MIAC, VLE *p*, and density, the same expression was introduced in the Table 10 footnote is used.

$$n_{a_1 c_1} = \frac{2\Delta^{a_1 c_1} n_{a_1} n_{c_1}}{V + \Delta^{a_1 c_1} n_t + \left[ V^2 + 2 \Delta^{a_1 c_1} V n_t + (\Delta^{a_1 c_1})^2 (n_{c_1} - n_{a_1})^2 \right]^{1/2}} \quad (17)$$

where $\Delta^{a_1 c_1}$ and $n_t$ are the equilibrium constant of an ion pair and the total number of ions, respectively. This expression (Eq (16)) is identical to that shown in Eq. (13), when only ion pairs can be formed. Hence, in this work, it is proposed to use the SAFT association term (Eq. (13)) for describing this contribution.

### 2.3.3. Computation of the association strength using the SAFT association term

Each pair of sites will require an association strength for which the original expression proposed by Wertheim theory is:

$$\Delta^{AB} = d_{ij}^3 g_{ij}^{hs} \kappa^{A_i B_j} \left[ \exp\left(\frac{\epsilon^{A_i B_j}}{kT}\right) - 1 \right] \quad (18)$$

where $\epsilon^{A_i B_j}$ is the energy of association, $\kappa^{A_i B_j}$ denotes the association volume and $g_{ij}^{hs}$ is the hard sphere *i-j* radial distribution function. In this work, we will keep this functionality for the solvent-solvent and the solvent-ion interactions but will reconsider it for ion-ion interactions. This association strength is in fact an apparent association constant expressed in volume units and should be related to the true association constant ($K_{AB}$). The originality of the association theory [70–72] is that the interactions are considered between sites of each molecule rather than between molecules. Therefore, it is considered that the formation of a pair ($A_i B_j$) is given by the interaction of site *A* of specie *i* ($A_i$) with site *B* of specie *j* ($B_j$). Thus, the equilibrium between sites can be expressed as:

$$A_i + B_j \rightleftharpoons A_i B_j \quad (19)$$

From Eq. (19), and considering cations site *A* and anions site *B*, the density-based association strength of the pair of sites ($\Delta^{A_i B_j}$) is given as:

$$\Delta^{A_i B_j} = \frac{\rho_{AB}}{\rho_A^+ \rho_B^-} \quad (20)$$

where $\rho_{AB}$, $\rho_A^+$ and $\rho_B^-$ are the density of pairs, density of cations and density of anions, respectively (the indices *i* and *j*, pointing to the molecular species that the site *A* and *B* belong to, are omitted for the sake of simplicity). On the other hand, the chemical theory considers a true association constant ($K_{AB}$) which, at equilibrium, is equal to the ratio of activities that are themselves related to the densities through $a_i = \rho_i / \rho_i^0 \gamma_i$, where $\rho^0$ is the unit density (1 mol/L).

$$K_{AB} = \frac{a_{AB}}{a_A a_B} = \frac{\rho_{AB} \gamma_{AB} \rho^0}{\rho_A^+ \gamma_A^+ \rho_B^- \gamma_B^-} = \frac{\rho_{AB} \rho^0}{\rho_A^+ \rho_B^-} \frac{\gamma_{AB}}{\gamma_A^+ \gamma_B^-} \quad (21)$$

The activity coefficients ($\gamma$) shown here are concentration (molarity) based. As can be seen from Eq. (21), the true association constant is expressed as a function of density and activity coefficients. As a result, the association strength and true association constants are related by the following equation:

$$\Delta^{AB} = \frac{K_{AB}}{\rho^0} \frac{\gamma_A^+ \gamma_B^-}{\gamma_{AB}} \quad (22)$$

which means that, in the equation of state framework, the $\Delta^{AB}$ (apparent association constant or association strength) can be obtained from the true association constant multiplied to a correction factor which is the
<figure><img src="image_5.png" /><figcaption>Fig. 2. Comparisons of the MIAC of (water + ethanol + NaCl) calculated using WS, BS, BT and BgT models at 298 K (lines are models result. Circles in black, gray, gold and green are experimental data related to pure water, 20 %, 40 %, and 80 % of ethanol, respectively).</figcaption></figure>

ratio of the concentration-based activity coefficients. It could be approximated by the following equation:

$$ \Delta^{AB} = \frac{K_{AB}}{\rho^0} \frac{\gamma_A^+ \gamma_B^-}{\gamma_{AB}} \approx K_{AB}(T) \beta_{AB} \quad (23) $$

where $K_{AB}(T)$ depends on the selected reference states but is generally considered as only a function of temperature and sometimes of the solvent composition, while $\beta_{AB}$ is a correction factor resulting from non-ideality. Eventually, Eq. (23) is rewritten follows to be used in the new framework:

$$ \Delta^{AB} = \Delta_{AB}^{0} \beta_{AB} \text{ where } \beta_{AB} = \frac{\gamma_{A}^{+} \gamma_{B}^{-}}{\gamma_{AB}} \quad (24) $$

where $\Delta_{AB}^0$ is the association constant of ions at infinite dilution ($x_{ion→0}$) in the solvent and $\beta_{AB}$ is the ratio of activity coefficients. The BiMSA theory provides a way to compute $\beta_{AB}$ which is also called the cavity function (see section (2.3.3.2)). Note that, for a given solvent, whereas the equilibrium constant of the chemical theory, $K_{AB}$, is only temperature-dependent, in the equation of state framework, both $\Delta_{AB}^0$ and $\beta_{AB}$ may be volume dependent (volume of solution). Consequently, even though $\Delta_{AB}^0$ is not explicitly salt concentration dependent, it can possibly vary with the salt content, simply because the molar volume decreases.

### 2.3.3.1. $\Delta_{AB}^0$ at infinite dilution.
Regarding the infinite dilution contribution, $\Delta_{AB}^0$ for ion pairing, two options are investigated. The first one is to use Eq. (18), where one can state that the correction factor $\beta_{AB}$ is equal to $g_{AB}^{hs}$, and what remains is $\Delta_{AB}^0$:

$$ \beta_{AB} = g_{AB}^{hs} \text{ and } \Delta_{AB}^{0} = d_{AB}^{3} \kappa^{AB} \left[ \exp\left(\frac{\epsilon^{AB}}{kT}\right) - 1 \right] \quad (25) $$

This is what was done in prior studies [33,34].

The second option is to use the Bjerrum theory [40]. This is an electrostatic theory that was introduced for the calculation of the association constant of ion pairs at infinite dilution. The Bjerrum association constant is obtained by solving the numerical integral expressed below:

$$ \Delta_{AB}^{0} = 4 \pi N_{A} (B \sigma_{AB}^{hs})^{3} \int_{2}^{b} (e^{x} x^{-4}) dx, \text{ with } B = \frac{2l_{B}}{\sigma_{ij}^{hs}} \quad (26) $$

$$ l_{B} = \frac{b_{c-a} |z_{a} z_{c}| e^{2}}{8 \pi \epsilon_{o} \epsilon_{r} k_{B} T} \quad (27) $$

where B is the reducing coefficient. $l_B$ is the Bjerrum length and $b_{AB}$ (later referred to as $b_{c-a}$ for cation-anion) is an adjustable parameter that is introduced for each ion pair (anion-cation) as a correction factor, allowing to differentiate the association strength with the type of ion pair. In this work, for the calculation of Bjerrum association strength, the
<figure><img src="image_6.png" /><figcaption>Fig. 3. Comparisons of the MIAC of (water + methanol + NaCl) calculated using WS, BS, BT and BgT models at 298 K (lines are models result. Circles in black, blue, gray, gold and green are experimental data related to pure water, 40 %, 60 %, 80 % and 100 % methanol, respectively). Data for 100 % methanol were taken from ref. [77] and, they were not included in regression (extrapolation results). The obtained AAD% for 100 % of Methanol using the WS, BS, BT and BgT models is 16.43, 3.24, 3.04 and 7.10, respectively.</figcaption></figure>

correlation introduced in previous work [34] is used:

$$ \Delta_{AB}^{0} = 4 \pi N_{A} (B\sigma_{AB}^{hs})^{3} e^{n_{1}|\ln(B-2)|^{3}+n_{2}|\ln(B-2)|^{2}+n_{3}\ln(B-2)+n_{4}(B-2)^{n_{6}}+n_{5}} \quad (28) $$

where the universal parameters are $n_1 = -0.03972$, $n_2 = -0.16739$, $n_3 = 0.6682$, $n_4 = 0.05678$, $n_5 = -1.16253$, $n_6 = 1.80043$. The parameters are obtained by a fit to the numerical integral in the range of $2.1 \le b \le 16$. The lower bound corresponds to $l_B = 1.05\sigma_{AB}^{hs}$, while the upper bound corresponds to $l_B = 8\sigma_{AB}^{hs}$. The upper bound defines an applicable range for $\epsilon_r$, which is meaningful for mixed-solvent electrolyte solutions, as $\epsilon_r$ of the co-solvent is usually much smaller than that of water.

### 2.3.3.2. Cavity function.
According to Eq. (29), to consider deviations from ideality in the association term (association of ions) a correction factor ($\beta_{AB}$) is used. It comes from the BiMSA model where it is called the cavity function:

$$ \Delta^{AB} = \Delta_{AB}^{0} \beta_{AB} \text{ where } \beta_{AB} = g_{AB}^{hs} g^{el} \quad (29) $$

where $\Delta_{AB}^{0}$ represents the association constant of ions at infinite dilution, and it can be obtained using the Bjerrum or Wertheim approach, experimental data or other approaches developed for the equilibrium

constant calculations of ion pairs. $g_{AB}^{hs}$ is the radial pair distribution function designating the short-range part of the interactions, in the primitive framework namely the contributions of hard sphere encountered for neutral species.

The value of $g_{AB}^{hs}$ starts from one (at infinite salt dilution) and reduces slightly with increasing salt concentration. However, in this work, it is assumed that this parameter is equal to one for the following reasons: first, to reduce the model complexity, and second, the impact of this expression was found to be negligible in the EoS framework. The second correction factor in Eq. (29) is $g^{el}$ which accounts for electrostatic ion screening. It is expressed as:

$$ g^{el} = \exp\left(\frac{2e^2(Z_aZ_c - Z_aZ_c')}{4\pi\epsilon_0\epsilon kT\sigma_p}\right) \quad (30) $$

This factor ($g^{el}$) considers the effect of long-range electrostatic interactions between pairing ions. The expressions for $Z$ can be found in Eq. (4).

## 2.4. List of modeling options

As a consequence of the above discussion, several models can be constructed and will be investigated in what follows. They are summarized in Table 3:
<figure><img src="image_7.png" /><figcaption>Fig. 4. Comparison of $\Delta_{AB}^0$ for sodium chloride salts in water using the WS and BS models in the temperature range of 298 K to 420 K.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Comparison of $\Delta_{AB}^0$ for sodium chloride salt in water-ethanol mixture using the WS and BS models at 298 K.</figcaption></figure>

Table 3 illustrates that in the ePPC-SAFT framework different models can be used for ion-ion association. For each of these ion-ion association models there are several choices for the RSP model. For each model, optimal parameters have been fitted as discussed in Section 2.5 and the results and discussion are presented in Sections 3 and 4.

## 2.5. Parameterization method

The solvent parameters of the PPC-SAFT model were obtained from previous works [33,73]. For parameterization of the ionic species and their interactions, two methods may be proposed: ion-specific and salt-specific. Under the salt-specific approach, the model parameters are exclusively valid for the corresponding salt. Conversely, in the ion-specific approach, the parameters apply to all salts featuring the

specific ion. The ion-specific method is favored as it simplifies extension to mixed-salt systems, which is in line with the final objective of this work. Furthermore, it should make it possible to yield physically meaningful parameter. However, obtaining these parameters requires regressing the parameters for a series of electrolyte solutions simultaneously. In this section the parameterization procedure for the chloride salts NaCl, KCl, RbCl and CsCl, in aqueous and mixed-solvent electrolyte solutions, is examined. It should be noted that, using the sBiMSA within the new framework does not add any adjustable parameter ($g^{el}$ is a predictive term) to the model. Therefore, the parametrization procedure and the number of adjustable parameters are the same as the ones used in previous works [33,34]. Data sourced from prior works on aqueous (Table 4) and mixed-solvent (Table 5) systems [33,38,41,62,63] are utilized in this work as well.

### 2.5.1. Adjustable parameters

In the same way as in previous work [34], the hard sphere ratio, $r^{HS}$ and the MSA ratio $r^{MSA}$ are defined as adjustable parameters for every ion, rather than the diameters themselves:

$$r^{HS} = \sigma_{HS} / \sigma_{Pauling} \qquad (31)$$

$$r^{MSA} = \sigma_{MSA} / \sigma_{HS} \qquad (32)$$

where the Pauling diameters are taken from ref [74]. The MSA diameters ($\sigma_{MSA}$) are those used in the BiMSA part of the model (Eq. (2)-6). The $\sigma_{HS}$ is the diameter used in the SAFT terms ($A^{hc}$ and $A^{disp}$). It is also used in the association term, as well as in the Wertheim and Bjerrum expressions (Eqs. (18) and (28)). The Born diameter, $\sigma_{Born}$, which is used in the Born term is not considered adjustable. It is taken from a previous work [61] where it was fitted on the solvation energy (see Table 7). For all the expressions that use the Wertheim equation (Eq. (18)), the numerical value of $\kappa^{AB}$ was systematically fixed to that of water, i.e. $\kappa^{AB} = 0.044394$. Also note that the dispersive interactions of the ions are entirely neglected, which means that their dispersive energy, $\epsilon$, is set to zero. Regarding the ion pairing parameters, we will distinguish between Wertheim ion-ion association, denoted as IP-Wertheim (Eq. (18)) and Bjerrum ion-ion association (Eq. (28)), denoted as IP-Bjerrum. In both cases, the solvation association energy (ion-solvent) is considered as an adjustable parameter. In practice, the CR1 combining rule [75] is used, such that the parameter to adjust is the “pseudo-unary” ion energy, $\epsilon_i^{AB}$:

$$\epsilon_{w-i}^{AB} = \frac{\epsilon_{w}^{AB} + \epsilon_{i}^{AB}}{2} \qquad (33)$$

where the subscript w points to water and the subscript i denotes an ion. When an alcohol is involved, then the ion-alcohol solvation strength is adjusted using a binary ion-alcohol parameter $w_{i-alcohol}$:

$$\epsilon_{i-alcohol}^{AB} = \frac{\epsilon_{i}^{AB} + \epsilon_{alcohol}^{AB}}{2} (1 - w_{i-alcohol}) \qquad (34)$$

In the IP-Wertheim case the ion-ion association energy is obtained from:

$$\epsilon_{c-a}^{AB} = \frac{\epsilon_{c}^{AB} + \epsilon_{a}^{AB}}{2} (1 - w_{c-a}) \qquad (35)$$

where the subscript c denotes a cation, and the subscript a an anion. In reality, the parameters are binary, yet it was found simpler to fix the cation and anion pseudo-unary parameters ($\epsilon_i^{AB}$) on the solvation (Eq. (33)), so that the binary parameter $w_{c-a}$ is used to adjust the ion pair association strength (Eq. (35)). Unary parameters refer to a property that is specific to the species, typically its diameter. On the other hand, binary parameters refer to an interaction. Note that all the binary parameters are symmetric, i.e. $w_{c-a} = w_{a-c}$, $b_{c-a} = b_{a-c}$. In the IP-Bjerrum case, the ion pairing association strength is computed from Eq. (28), and the binary parameter $b_{c-a}$ (Eq. (27)) is used to adjust the strength. It
<figure><img src="image_9.png" /><figcaption>Fig. 6. 6a shows the molar volume trend in equation of state framework using BS and BT models. In EoS framework volume (molar volume) has always negative slope using any type of model. 6b and 6c compare the dielectric constant (RSP) and $\Delta_{AB}^0$ of sodium chloride salt in water using the BT and BS models at 298 K.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. 7a shows $g^e$ shape using BgT model for sodium chloride salts in water at 298 K. 7b and 7c present the obtained association constant of ions for sodium chloride salt in water at 298 K using the BgS and BgT models respectively.</figcaption></figure>

must be noted that, according to Table 6, for WS model the Wertheim model (IP-Wertheim) is used for ion-ion association, and therefore $w_{c-a}$ parameters, while for the rest (BS, BT, BgS and BgT) the Bjerrum model (IP-Bjerrum) and therefore the $b_{c-a}$ parameters.

## 2.5.2. Parameterization of aqueous chloride salts solutions

In a first stage the chlorides are considered in an aqueous solution. Therefore, according to Table 6, for chloride salts in water 15 parameters are regressed simultaneously. Note that the number of adjustable parameters is equal for the five models. The only difference relates to the type of ion-ion interaction parameter in which is ($b_{a-c}$) for Bjerrum and ($w_{c-a}$) for Wertheim. As in our previous works, it was chosen to keep the same value of $r^{HS}$ for all ions (to keep the anion and cation parameter in the same orders). The following objective function is used for the model parametrization of aqueous chloride salts solutions:

$$ \text{OF} = \frac{\sum_{i=1}^{n_{dp}^{\pm}} \left( \frac{\gamma_{i,cal}^{\pm} - \gamma_{i,exp}^{\pm}}{\gamma_{i,exp}^{\pm}} \right)^2}{n_{dp}^{\gamma^{\pm}}} + \frac{\sum_{i=1}^{n_{dp}^{\rho}} \left( \frac{p_{i,cal} - p_{i,exp}}{p_{i,exp}} \right)^2}{n_{dp}^{\rho}} + \frac{\sum_{i=1}^{n_{dp}^{\rho}} \left( \frac{\rho_{i,cal} - \rho_{i,exp}}{\rho_{i,exp}} \right)^2}{n_{dp}^{\rho}} \quad (36) $$

where $n_{dp}$ is the number of data points for the property, $\gamma^{\pm}$ is the MIAC, $\rho$ is the vapour pressure and $\rho$ is the molar density. All parameters of the four models are shown in Table 7 (unary and pseudo unary) and Table 8 (binary). In practice, the BgS approach never delivered acceptable results (as will be further detailed in Section 4), which is why it is skipped in the first section. Table 7 shows in fact for WS, BT models, trends are

observed in the ion-ion ($b_{a-c}$, $w_{c-a}$) and cation-water ($\epsilon_{w-i}^{AB}$) association. When considering cations of increasing size, a consistent decrease in the ion-water association strength and a consistent increase in the cation-anion association strength are observed for WS and BT models. Regarding the $r^{MSA}$, a minimum is often observed for this ratio for the $Rb^+$ cation ($K^+$ in the case of BgT). The MSA diameter is generally considered to be the solvated diameter, and it may therefore be larger for smaller cations that are more strongly solvated. The deviations are discussed in Section 3 Table 8.

## 2.5.3. Parameterization of mixed-solvent chloride salts solutions

In a second stage, the parameterization is performed on mixed solvents data (water + methanol and water + ethanol). It was chosen not to work with binary alcohol-salt data, because these data only cover a very small concentration range. Hence, they do not contain sufficient information.

The ion-specific parametrization procedure for mixed-solvent chloride salts solutions has been introduced and detailed in previous work as well [34]. Here the same procedure is followed. For the mixed solvent systems only the association energy interaction parameter of ions with alcohol ($w_{ion-alcohol}$) as shown in Eq. (34) is regressed and the rest of the parameters are taken from the aqueous system. The objective function for the mixed-solvent system is:
<figure><img src="image_11.png" /><figcaption>Fig. 8. Comparison of ion-solvent and ion-ion association strength at infinite dilution (Δ₀ₐₗ) using the four models for chloride salts in water at 298 K (ion-ion association in blue and ion-solvent association in black). Water-water association constant obtained from all models is the same (water-water = 1.66 10⁻³ m³/mol). For water-Cl⁻, the obtained values for association constant using the WS, BS, BT and BgT models are 2.30 10⁻³, 1.65 10⁻², 1.73 10⁻² and 3.84 10⁻³ m³/mol respectively.</figcaption></figure>

$$ \text{OF} = \frac{\sum_{l=1}^{n_{\text{dp}}^{\pm}} \left( \frac{y^{\pm}_{i, \text{cal}} - y^{\pm}_{i, \text{exp}}}{\gamma^{\pm}_{i, \text{exp}}} \right)^2}{n_{\text{dp}}^{\pm}} + \frac{\sum_{l=1}^{n_{\text{dp}}^{\text{p}}} \left( \frac{p_{i, \text{cal}} - p_{i, \text{exp}}}{p_{i, \text{exp}}} \right)^2}{n_{\text{dp}}^{\text{p}}} + \frac{\sum_{l=1}^{n_{\text{dp}}^{\text{y}}} \left( y_{i, \text{cal}} - y_{i, \text{exp}} \right)^2}{n_{\text{dp}}^{\text{y}}} + \frac{\sum_{l=1}^{n_{\text{dp}}^{\text{y}}} \left( \frac{p_{i, \text{cal}} - p_{i, \text{exp}}}{p_{i, \text{exp}}} \right)^2}{n_{\text{dp}}^{\text{y}}} \quad (37) $$

where y is the mole fraction of the alcohol in the vapour phase. The γ± point here to mean ionic activity coefficients (MIAC), that are molality based and use a mixed solvent reference state. The alcohol-ion binary interaction parameters (BIPs) for the four models and the two mixed solvent systems are presented in Table 9.

### 3. Comparison of the models

In this section results from WS, BS, BT and BgT models for chloride salts in aqueous and mixed solvents systems are compared and discussed.

#### 3.1. Chloride salts in aqueous system

Table 10 shows the average absolute percentage deviations (AADs) and maximum absolute percentage deviations (MADs) for MIAC, vapour

pressure and density of WS, BS, BT and BgT models for aqueous chloride salts solutions. The use of the maximum deviation along with the average deviation is of great interest in order to show whether the model follows the expected trend. Focusing on MIAC, the most precise model is clearly the WS model, both on average and regarding maximum deviations. For vapour pressures, the BT model is slightly more accurate, and it is also acceptable for MIAC. For density, BgT is probably best, although all models provide on average good results. Globally speaking, the least accurate is the BS model. It appears that the combination of the Bjerrum ion pairing with the Schreckenberg RSP makes it difficult to fit the experimental data, as is further discussed below.

#### 3.2. Chloride salts in mixed solvent systems

Table 11 shows the average absolute percentage deviations (AADs) and maximum absolute percentage deviations (MADs) for MIAC, bubble pressure (p), vapour composition (y), and density of the four models WS, BS, BT and BgT for chloride salts in mixed solvent systems (water-ethanol and water-methanol). The deviations are much larger than those observed with aqueous systems, which illustrates the complexity of the mixed solvent electrolyte modeling. For densities, average deviations often remain close to 1 %, but are higher for the methanol solvent. All models perform very similarly. The same can be stated for the vapour-liquid equilibrium properties (bubble pressure and vapour composition).
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison of association constant of ions ($\Delta_{AB}^0$) using the WS, BS, BT and BgT models for chloride salts in water ethanol mixture at 298 K. W-W, W-E and E-E represent the water-water, water-ethanol and ethanol-ethanol association strengths respectively. Cl-E and Na-E represent ethanol-ions association strengths. Cl-W and Na-W represent water-ions association strengths. Na-Cl represents the ion-ion association strength.</figcaption></figure>

The MIAC, which is often considered as the most difficult property to describe accurately, indeed shows the highest maximum deviations. Deviations of the order of 5 % on average are considered acceptable, with maximum deviations above 10 %. Other authors found in fact similar deviations [11,30,32–34,62]. The WS and BgT models yield the highest deviation with respect to the experimental data. For water-methanol system, the BT model and for water-ethanol system, the BS model are the most precise models in terms of accuracy.

Regarding the MIAC, the results for NaCl in ethanol-water and methanol-water mixtures are shown in Figs. 2 and 3. For these systems, a large number of data are available in an extensive range (Table 5). All models are very accurate for pure water, but deteriorate with increasing alcohol content, more so with ethanol than with methanol. Note that no binary salt-alcohol data have been used during fitting. The figures show such data (100 % alcohol) and the quality of the predictions. As can be seen, using the BT model the best result (AAD% = 3.04 is obtained which shows the capability of this model for extrapolation to higher methanol concentration. Note also that almost all parameters are fitted on pure water solvent. Only one parameter, the ion-alcohol solvation energy, is fitted on the mixed solvent. Yet, the results show that this is not sufficient.

# 4. Results and discussion

The five model combinations (WS, BS, BT, BgS and BgT) are discussed in this section. The discussion focuses essentially on the model capacity to describe ion solvation and ion pair formation using the association and MAL terms, respectively. The first section focuses on ion pairing; next, we compare the different association strengths and finally discuss the contributions to the activity coefficient calculations.

## 4.1. Ion pairing association strength

In this section a comprehensive discussion of the ion-ion association models (Bjerrum or Wertheim) and the effect of RSP on ion-ion association strength is made.

### 4.1.1. Comparison of Bjerrum and Wertheim for the calculation of $\Delta_{AB}^0$

The choice between the Bjerrum or the Wertheim theory to describe ion pairing is the first question discussed in this work. Previous investigations [34] have pointed to the fact that the use of ion pairing with the Wertheim functionality (WS) provided much better results than that of Bjerrum for aqueous systems. Yet, it was unable to be extended to mixed solvents unless “effective” RSP values were introduced. Here, a deeper analysis of the WS and BS models is proposed using chloride salts in aqueous and mixed solvents systems. Both models use the same RSP
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Contributions to lnγ± of the terms for chloride salts in water using the WS, BS, BT and BgT models at 0.1 MPa and 298 K. The sum is shown in dashed blue, and the color code for the contributions is as follows; Hard chain (HC) in blue, dispersion (Disp) in orange, Association (Assoc) in gray, Polar in gold, Born in green, MSA in yellow and MAL in dark blue.</figcaption></figure>

model so that comparing these two models effectively allows comparing
the IP-Wertheim and IP-Bjerrum approaches.

4.1.1.1. *Chloride salts in aqueous systems.* In Fig. 4, the Δ₀<sup>AB</sup> values for sodium chloride in water are plotted using the WS and BS models. The figure reveals that, according to the Wertheim theory (WS model), as

temperature increases, Δ₀<sup>AB</sup> values decrease. Conversely, according to the Bjerrum theory (BS model), Δ₀<sup>AB</sup> increases with temperature. Yet, it is known that the ion pair concentration increases by raising the temperature [78–80], which is indeed what the BS model shows in contrast to WS. This is the result of the fact that the Wertheim association (Eq. (18)) is independent of the dielectric constant (RSP). However, for lower
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Contributions to ln(γ<sub>±</sub>) of the terms for sodium chloride in pure ethanol using the WS, BS, BT and BgT models at 0.1 MPa and 298 K. The sum is shown in dashed blue, and the color code for the contributions is as follows; Hard chain (HC) in blue, Dispersion (Disp) in orange, Association (Assoc) in gray, Polar in gold, Born in green, MSA in yellow and MAL in dark blue.</figcaption></figure>

temperatures (close to 298 K, where the parameters have been fitted),
both association models are very close.

4.1.1.2. Chloride salts in mixed solvents systems. In Fig. 5, the Δ₀<sup>AB</sup> values for sodium chloride in mixed solvent system (ethanol-water) are plotted using WS and BS models. According to this figure, Δ₀<sup>AB</sup> slightly decreases with an increase in ethanol concentration when using the Wertheim theory (WS model). Conversely, according to the Bjerrum theory (BS model), Δ₀<sup>AB</sup> strongly increases with increasing ethanol concentration. It is well-established that the presence of alcohols in the system signifi- cantly reduces the dielectric constant, leading to a notable increase in ion pair concentration [49,50]. Hence, the BS model clearly behaves more consistently with the physical expectation, which is related to the fact that BS incorporates the impact of the dielectric constant.

As a result, it can be concluded that the Wertheim theory is not a
good solution for the calculation of ion-ion association for aqueous and
mixed solvent systems specially when the concentration of ion pairs is
significant (high temperature or alcohol concentration).

4.1.2. Impact of RSP on Δ₀ᵃᵇ

In contrast to the Wertheim equation for the calculation of Δ₀<sup>AB</sup>, the Bjerrum theory is sensitive to the RSP. Hence, the choice of the functionality is essential. Its impact is now investigated by comparing the BT and BS models.

At first, it must be recognized that the molar volume of the salt solute decreases upon addition of salt (Fig. 6a). Using the Schreckenberg model, the RSP is explicitly volume-dependent, such that this volume decrease leads to a reduction in the RSP (as depicted in Fig. 6b). Note that this is indeed the experimentally expected behavior. Yet, consequently, because the Bjerrum association strength depends on this RSP, the value of Δ₀<sup>AB</sup> in the BS model increases (Fig. 6c). However, based on Eqs. (23) or (24), Δ₀<sup>AB</sup> represents the association constant of ions at infinite dilution and should therefore not be affected by the ion concentration. According to Eq. (29), only g<sup>el</sup> drives the association constant of ions at high concentration. The BS model thus contains an inconsistency.

In order to have a salt concentration independent Δ₀ᵃᵇ for the Bjerrum model, the RSP must be independent of both volume and ion concentration (only a function of temperature according to Eq. (23)). This is why we introduced the Fürst-Zuo (Eq. (10)) and Maryott-Smith (Eq. (11)) RSP approach in the BT model. This leads to a RSP that is independent of salt concentration (Fig. 6b), which is more consistent with the primitive approach that is constructed on a framework that considers the solvent through a constant dielectric constant.

4.1.3. Impact of the cavity function (g^el) on Δ^0_{AB}

The impact of salt concentration on the association strength is
described using g^el according to Eq. (29). To investigate the effect of RSP
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Contributions to difference between lnγ<sub>alc</sub> (logarithmic alcohol activity coefficient) at ion mole fraction 0.03 and in the salt-free solvent, for (water + ethanol + NaCl) at 298.15 K and 0.1 MPa. This difference is positive (see dashed line for the sum). The sum is shown in dashed blue (lnγ<sub>alc</sub>), and the color code for the contributions is as follows; Hard chain (HC) in blue, Dispersion (Disp) in orange, Association (Assoc) in gray, Polar in gold, Born in green, MSA in yellow and MAL in dark blue. The pictures show the four models WS, BS, BT and BgT.</figcaption></figure>

on g<sup>el</sup>, BgS and BgT models are compared. It should be noted that the
results using the BgS model are not acceptable, not only in terms of
accuracy but also in terms of parameter consistency (very large r<sup>MSA</sup>
values are obtained from the regression), and they are presented in
supplementary information. The shape of g<sup>el</sup> with ion concentration in
water for the BgT model is shown in Fig. 7a: it is a decreasing function of
the salt content, and its value at infinite dilution is one. It decreases very
quickly at low salt concentration and stabilizes at a value close to 0.4
m<sup>3</sup>/mol. Consequently, the actual association strength (based on Eq.
(29)), pictured in Figs. 7a and b, has either a monotonic decreasing
behavior (BgT, 7c) or a minimum (BgS, 7b). For the BgT model (Fig. 7c),
using the Zuo & Fürst RSP, the association constant Δ has the shape that
is expected using the BiMSA approach. In contrast, the shape of the
association strength using the BgS model (Fig. 7b) shows a minimum.
This could be the reason why the BgS model fails in our EoS framework.

## 4.2. Comparison and impact of the association strengths: solvation and ion pairing

Collins [81] points to the importance of describing solvation more
explicitly than is usually done in electrolyte models. Wexler et al. [42,
82–84] published a series of papers stating that the non-idealities in
electrolyte solutions can essentially be described using pseudo-chemical
solvation theories. In this work, the pseudo-chemical theory is taken into
account through the association contributions. As pointed out above,
three types of association interactions can be distinguished:
solvent-solvent, solvent-ion (that we call solvation) and ion-ion (that we
call ion pairing). Comparing the magnitudes of the association strengths

at zero salt concentration (Δ<sup>0</sup><sub>AB</sub>), we can evaluate qualitatively the value
of the different models.

Fig. 8 compares the cation-solvent and ion-ion association strengths at zero salt content, obtained using the WS, BS, BT and BgT models for chloride salts in water. The black lines show the solvation strength, and the blue are the ion pairing strengths. As can be seen in this figure for WS and to a lesser extend for BT (except for Na<sup>+</sup>), a trend is observed for hydration and ion pairing strengths. From left to right (NaCl to CsCl), the hydration reduces and conversely, the ion pairing increases. This is consistent with the observation by Vaque Aura et al. [1] who noticed that in the case of chloride salts there is a consistent decrease in MIAC from NaCl to CsCl, which can be related to the law of matching water affinities developed by Collins [85]. According to this theory, when considering the series of chloride salts from LiCl to CsCl, the solvation strength should decrease, and the ion pairing should increase.

The obtained results using the WS, BT and BgT models also show that anion-water association strength is greater than that of cation-water. Anion-water association strength is in the range of 10<sup>-5</sup> to 10<sup>-4</sup> (m<sup>3</sup>/mol) and for cation-water, it is in the range of 10<sup>-3</sup> to 10<sup>-2</sup> (m<sup>3</sup>/mol). Collins [81] also explained why this may be consistent: Firstly, quantum mechanical calculations indicate that the anions which interact with the hydrogen atom of water, allow intra-shell hydrogen bonding of the solvating waters, whereas cations, which interact with the oxygen atom of water, do not [86]. Second, charge transfer between ion and solvent characterizes strong hydration, and water accepts negative charge from anions more efficiently than it donates negative charge to cations.

In order to illustrate the importance of ion pairing in mixed solvent
solutions, Fig. 9 shows how the infinite dilution association strengths
depend on the solvent composition for the case of mixtures with ethanol.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Contributions to the difference between lnγ<sub>w</sub> (logarithmic water activity coefficient) at ion mole fraction of 0.03 and its salt-free solvent, for the (water + ethanol + NaCl) at 298.15 K and 0.1 MPa. This difference is negative (see dashed line for the sum). The sum is shown in dashed blue (lnγ<sub>w</sub>), and the color code for the contributions is as follows; Hard chain (HC) in blue, Dispersion (Disp) in orange, Association (Assoc) in gray, Polar in gold, Born in green and MSA in yellow and MAL in dark blue. The pictures show the four models WS, BS, BT and BgT.</figcaption></figure>

The solvent-solvent hydrogen bonding strengths are almost identical in all cases (water-water, ethanol-ethanol and water-ethanol). They only depend on the solution molar volume which varies with the solvent composition. Similarly, the ion-solvent association strengths almost do not change with solvent composition. The main observation of these plots is the drastic change in ion-ion association strength: in pure water solvent (on the left), there is almost no ion pairing, while in pure ethanol (on the right), the ion pairing strength is much larger for the Bjerrum-based models. The Wertheim-based model is incapable of describing this behavior. On a logarithmic representation, all Bjerrum-based models yield very similar results: the ethanol-chloride solvation strength is very similar to the ethanol-ethanol hydrogen bonding strength, and the Na⁺-ethanol solvation is much smaller. The Na⁺-water solvation strength is very similar to that of Na⁺-ethanol, but the Cl⁻-water strength is much stronger than that of Cl-ethanol, which is expected due to the high polarity of water.

## 4.3. Contributions of different terms

In this section, the contributions of the terms to the MIAC and solvent activity coefficients in the aqueous and mixed-solvent electrolyte solutions are detailed and discussed.

### 4.3.1. Contributions to MIAC

The mean ionic activity coefficient (γ<sub>±</sub>) is defined as the ratio of fugacities in the mixture and in the reference state (*), here the ion free mixed solvent. Since the equation of state provides mole fraction-based activity coefficient, we will analyze this value, which is easy to trans-

form into molality base:

$$RT\ln\gamma_{\pm} = \sum \left( \frac{\partial A^x}{\partial n_{ion}} - \frac{\partial A^{x+}}{\partial n_{ion}} \right) - RT\ln \frac{Z}{Z^*} \quad (38)$$

where X refer to the different contributions from the EoS introduced in Eq. (1).

Fig. 10 shows the contributions to the mole fraction based mean ionic activity coefficient, lnγ<sub>±</sub> of terms for chloride salts in water (when we have minimum ion pairing) using the WS, BS, BT and BgT models. According to this figure in all models, the association term (assoc) is the dominant contribution and is positive. It describes how the structure of the solution changes upon addition of salt, and therefore replaces the so-called “Born” term that describes only the electrostatics of the ion-solvent interactions. This positive term is counterbalanced by the electrostatic ion-ion interaction, provided here by MSA. This is the expected behavior. Among the other contributions, it is worth pointing to the large negative impact of the hard sphere (HC) term, especially for the large cation (Cs⁺). The other contributions (including Born) are small.

Notice that in previous works [33,34], association contained all types of interactions, while here, the solvent-ion contribution is distinct from ion pairing, called MAL. It shows the important impact of structural solvation on the ionic activity coefficients. The MAL term is always negative and negligible in comparison to the association term, as expected in aqueous solvents. Notice that the Born term contribution is exactly zero when the RSP model is not salt-concentration dependent (BT and BgT) and very small otherwise.

Fig. 11 shows the contributions to the mole fraction based lnγ<sub>±</sub> of the EoS terms for sodium chloride in pure ethanol (when we have
maximum ion pairing) using the WS, BS, BT and BgT models. Here, the different models perform very differently: as already pointed out, the WS model yields a very small ion pairing strength, which is here visible through the negligible MAL term. As a consequence, it is impossible for this model to describe the experimental trend, since the negative MSA term has no parameter that can be fitted on alcohol containing solvents. Only a drastic reduction of the solvation strength can help decrease the total value of MIAC, which is also seen in Fig. 9 (WS has the smallest ion-alcohol interaction strength). In contrast, for Bjerrum based models (BS, BT and BgT), the MAL term is now strongly negative in absolute value and contributes significantly to the decrease of the ionic activity co-efficients. Unlike for the aqueous system, for the mixed solvent systems MAL is the dominant contribution, together with MSA that also strongly decreases as a result of the lower RSP of the alcohol. The solvation (association) remains a positive term, but much smaller than it was for aqueous systems. For the BgT model, even at high alcohol content, MAL does not increase as much as the other Bjerrum based models. This is because in this model, the g<sup>el</sup> term reduces the association constant of ions at high salt or alcohol concentrations.

### 4.3.2. Contributions to solvent activity coefficient

The activity coefficient of the solvents in salt-containing solutions is a measure of the salting-in or salting-out effect: when the activity co-efficient rises as a consequence of salt addition, the solvent becomes more volatile, which points to salting out; in contrast, when the salt activity coefficient drops, salting-in occurs. In the case of water-alcohol systems, water is salted in, because of what is generally called a reduction in “free water”. In contrast, the alcohol is then salted out because the solvation of alcohol is lower than that of water.

Figs. 12 and 13 then show the impact of the salt on the solvent activity coefficients. It is constructed using the same equation as above:

$$RT\Delta ln\gamma_i = \sum \left( \frac{\partial A^x}{\partial n_i} - \frac{\partial A^{x-x}}{\partial n_i} \right) - RT\ln \frac{Z}{Z^*} \quad (39)$$

where X refer to the different contributions from the EoS introduced in Eq. (1) and the asterisk point to the salt-free mixed solvent solution. The index i refers to the solvent, either alcohol or water. The plots show the difference between the derivatives for a salt-containing mixed solvent with 3 % mole fraction and the salt-free mixtures.

The models that are developed in this work make it possible to analyze the contribution of each term to this property (Figs. 12 and 13). Interestingly, the positive sign of the alcohol activity coefficient variation when adding salts is essentially driven in the positive direction by the dispersion term. This can be explained by the fact that in this work, the ions do not interact at all through dispersion. Hence, when salt is added, the dispersion energy (an attractive energy) decreases, resulting in a repelling effect for the solvent molecules. It is balanced by ion solvation and hard sphere terms that are negative, but globally, the presence of ions results in an increase of the alcohol activity coefficient. The picture is the same for both BT and BgT. Using the Schreckenberg RSP, the hard sphere contribution becomes positive, and the association term is still more negative.

Fig. 13, that illustrates the salting in of water (the sum is negative) shows also a positive contribution of the dispersive term, but this time, we observe a very strong negative contribution of both the Born and the association terms. These terms point to the solvation: Born describes its electrostatic part, and the association describes its structuring part (water is more strongly bound as a consequence of the presence of ions), and this is even more so in a solution that contains more alcohol.

## 5. Conclusion

In this work, the Binding-MSA model was implemented within the ePPC-SAFT framework, and the new model (BiMSA-ePPC-SAFT) was

applied to chloride salts in both aqueous and mixed-solvent systems different association strengths models were investigated for ion pairing.

Regarding the debates and concerns that were introduced at the beginning of this article, a few comments and conclusions are made:

* Which model (Wertheim or Bjerrum) best describes ion pairing in solution?

The Wertheim theory is unable to describe the trend of ion pairing with either temperature or solvent composition. In contrast, the Bjerrum based models are able to correlate with the experimental data of both aqueous and mixed solvent systems (specifically mixed solvent systems except BgT model), because of its dependence on dielectric constant. As a result, it could be concluded that the Bjerrum approach is more reliable or is able to describe the ion pairing phenomenon much better than the Wertheim theory.

* What is the best RSP functionality?

One important conclusion of this work was obtained from the comparison of different types of RSP models (in particular its volume-dependence), and more specifically their effect on the ion pairing behavior. It appears that, using the Bjerrum approach for ion pairing strength, the best RSP functionality is the one that does not decrease with salinity (neither explicitly, nor implicitly through its volume dependence).

* What type of interactions (ion-ion, ion-solvent, solvent-solvent, solvent-cosolvent or co-solvent-cosolvent) are more important for the modeling of aqueous and mixed solvent systems?

It was observed that, for an aqueous system, ion-solvent and solvent-solvent associations have the strongest impact on the mean ionic activity coefficient. Ion-ion interactions are taken into account through an electrostatic term (MSA), while solvation (ion-solvent interactions) is more impacted by structural effects than by electrostatic ones (as also pointed out by Collins [81,85]). However, in non-aqueous solvent systems, ion pairs are formed, which needs to be taken into account in the model through an additional MAL term.

Overall, the BiMSA model was successfully implemented within the ePPC-SAFT framework. Using the BS and BT models a systematic improvement was found for mixed solvent systems in comparison to WS model. The BgT model was not as accurate as anticipated in the case of mixed solvents. In other words, adding the cavity function (g<sup>el</sup>) to consider the non-ideality of ion-ion association in high salt concentrations did not improve the accuracy of the model. We believe that this point will deserve further investigation in the future. It was also shown that the model makes it possible to visualize the relative impact of the various physical interactions on activity coefficients, including the salting-in or salting-out phenomenon. It is also our intention to extend this approach in future investigations so that other salts and other types of phase equilibria are included.

### CRediT authorship contribution statement

Abtin Raeispour Shirazi: Writing – review & editing, Writing – original draft, Software, Methodology, Investigation. Fufang Yang: Validation, Software, Methodology, Data curation. Tri Dat Ngo: Visualization, Validation, Software, Methodology, Data curation. Nicolas Ferrando: Writing – review & editing, Validation, Methodology. Olivier Bernard: Writing – review & editing, Validation, Supervision, Methodology, Investigation. Jean-Pierre Simonin: Writing – review & editing, Validation, Supervision, Project administration, Methodology, Investigation. Jean-Charles de Hemptinne: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Conceptualization.
Declaration of competing interest

The authors declare the following financial interests/personal re-
lationships which may be considered as potential competing interests:
Jean-Charles de Hemptinne reports financial support was provided by
IFP Energies nouvelles. If there are other authors, they declare that they
have no known competing financial interests or personal relationships
that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgement

Financial support from electrolyte thermodynamic (Elether) chair is
gratefully acknowledged.

Supplementary materials

Supplementary material associated with this article can be found, in
the online version, at doi:10.1016/j.fluid.2024.114231.

References

[1] S. Vaque Aura, J.-S. Roa Pinto, N. Ferrando, J.-C.de Hemptinne, A.ten Kate, S. Kuitunen, N. Diamantonis, T. Gerlach, M. Heilig, G. Becker, M Brehelin, Data Analysis for Electrolyte Systems: A Method Illustrated on Alkali Halides in Water, J. Chem. Eng. Data 66 (8) (2021) 2976-2990, https://doi.org/10.1021/acs.jced.1c00105.

[2] Z.-G. Li, Y.-Q. Sun, W.-L. Zheng, H. Teng, Z.-L. Xiu, A novel and environment-
friendly bioprocess of 1,3-propanediol fermentation integrated with aqueous two-
phase extraction by ethanol/sodium carbonate system, Biochem. Eng. J. 80 (2013)
68-75, https://doi.org/10.1016/j.bej.2013.09.014.

[3] A. Jouyban, Handbook of Solubility Data For Pharmaceuticals, CRC Press, United States of America, 2010.

[4] O. Lefebvre, R. Moletta, Treatment of organic pollution in industrial saline wastewater: a literature review, Water Res. 40 (20) (2006) 3671-3682, https://doi.org/10.1016/j.watres.2006.08.027.

[5] P. Zhang, T. Yokoyama, O. Itabashi, T.M. Suzuki, K. Inoue, Hydrometallurgical process for recovery of metal values from spent lithium-ion secondary batteries, Hydrometallurgy 47 (2-3) (1998) 259-271, https://doi.org/10.1016/S0304-386X (97)00050-9.

[6] S. Naseri Boroujeni, X. Liang, B. Maribo-Mogensen, G.M Kontogeorgis, Comparison of Models for the Prediction of the Electrical Conductivity of Electrolyte Solutions, Ind. Eng. Chem. Res. 61 (8) (2022) 3168-3185, https://doi.org/10.1021/acs.iecr.04365.

[7] J.R. Loehe, M.D. Donohue, Recent advances in modeling thermodynamic properties of aqueous strong electrolyte systems, AIChE J. 43 (1) (1997) 180-195, https://doi.org/10.1002/aic.690430121.

[8] M.L. Pinsky, K. Takano, Property Estimation for Electrolyte Systems. Computer Aided Property Estimation for Process and Product Design, Elsevier, 2004, pp. 181-204.

[9] G.M. Kontogeorgis, B. Maribo-Mogensen, K. Thomsen, The Debye-Hückel theory and its importance in modeling electrolyte solutions, Fluid. Phase Equilib. 462 (2018) 130-152, https://doi.org/10.1016/j.fluid.2018.01.004.

[10] G.M. Kontogeorgis, A. Schlaikjer, M.D. Olsen, B. Maribo-Mogensen, K. Thomsen, Solms N. von, X Liang, A Review of Electrolyte Equations of State with Emphasis on Those Based on Cubic and Cubic-Plus-Association (CPA) Models, Int. J. Thermophysics. 43 (2022) 4, https://doi.org/10.1007/s10765-022-02976-4.

[11] N. Novak, G.M. Kontogeorgis, M. Castier, I.G. Economou, Mixed solvent electrolyte solutions: a review and calculations with the eSAFT-VR Mie equation of state, Ind. Eng. Chem. Res. 62 (34) (2023) 13646-13665, https://doi.org/10.1021/acs.iecr.3c00717.

[12] S. Ahmed, Thèse de doctorat, Sorbonne Université, 2019.

[13] C. Held, Thermodynamic gE models and equations of state for electrolytes in a water-poor medium: a review, J. Chem. Eng. Data 65 (11) (2020) 5073-5082, https://doi.org/10.1021/acs.jced.0c00812.

[14] H. Zhao, M.C. dos Ramos, C. McCabe, Development of an equation of state for electrolyte solutions by combining the statistical associating fluid theory and the mean spherical approximation for the nonprimitive model, J. Chem. Phys. 126 (24) (2007) 244503, https://doi.org/10.1063/1.2733673.

[15] G. Das, S. Hlushak, M.C. dos Ramos, C. McCabe, Predicting the thermodynamic properties and dielectric behavior of electrolyte solutions using the SAFT-VR+DE equation of state, AIChE J. 61 (9) (2015) 3053-3072, https://doi.org/10.1002/aic.14909.

[16] S. Herzog, J. Gross, W. Arlt, Equation of state for aqueous electrolyte systems based on the semirestricted non-primitive mean spherical approximation, Fluid. Phase Equilib. 297 (1) (2010) 23-33, https://doi.org/10.1016/j.fluid.2010.05.024.

[17] J.L. Lebowitz, J.K. Percus, Mean spherical model for lattice gases with extended hard cores and continuum fluids, Phys. Rev. 144 (1) (1966) 251-258, https://doi.org/10.1103/PhysRev.144.251.

[18] P. Debye, E. H. I. abaissement du point de congelation et phenomenes associes, Phys. Z. 24 (1923) 185-206. French.

[19] M. Born, Volumen und Hydratationswärme der Ionen, Zeitschrift für Physik 1 (1) (1920) 45-48, https://doi.org/10.1007/BF01881023.

[20] K.E. Gubbins, C.H. Twu, Thermodynamics of polyatomic fluid mixtures—I theory, Chem. Eng. Sci. 33 (7) (1978) 863-878, https://doi.org/10.1016/0009-2509(78)85176-8.

[21] D. NguyenHuynh, J.-P. Passarello, P. Tobaly, J.-C.de Hemptinne, Application of GC-SAFT EOS to polar systems using a segment approach, Fluid. Phase Equilib. 264 (1-2) (2008) 62-75, https://doi.org/10.1016/j.fluid.2007.10.019.

[22] P.K. Jog, S.G. Sauer, J. Blaesing, W.G. Chapman, P.K. Jog, S.G. Sauer, J. Blaesing, W.G. Chapman, Application of dipolar chain theory to the phase behavior of polar fluids and mixtures, Ind. Eng. Chem. Res. 40 (2001) 4641-4648.

[23] B. Maribo-Mogensen, G.M. Kontogeorgis, K. Thomsen, Comparison of the Debye-Hückel and the Mean Spherical Approximation Theories for Electrolyte Solutions, Ind. Eng. Chem. Res. 51 (14) (2012) 5353-5363, https://doi.org/10.1021/ie2029943.

[24] J.-P. Simonin, O. Bernard, Insight into the ionic atmosphere effect: comparison of theories for electrolytes at the primitive level, Fluid. Phase Equilib. 571 (2023) 113805, https://doi.org/10.1016/j.fluid.2023.113805.

[25] B. Maribo-Mogensen, K. Thomsen, G.M. Kontogeorgis, An electrolyte CPA equation of state for mixed solvent electrolytes, AIChE J. 61 (9) (2015) 2933-2950, https://doi.org/10.1002/aic.14829.

[26] S. Ahmed, N. Ferrando, J.-C.de Hemptinne, J.-P. Simonin, O. Bernard, O Baudouin, Modeling of mixed-solvent electrolyte systems, Fluid. Phase Equilib. 459 (2018) 138-157, https://doi.org/10.1016/j.fluid.2017.12.002.

[27] G. Das, M.C. dos Ramos, C. McCabe, Predicting the thermodynamic properties of experimental mixed-solvent electrolyte systems using the eSAFT-VR+DE equation of state, Fluid. Phase Equilib. 460 (2018) 105-118, https://doi.org/10.1016/j.fluid.2017.11.017.

[28] J.M. Schreckenberg, S. Dufal, A.J. Haslam, C.S. Adjiman, G. Jackson, A. Galindo, Modelling of the thermodynamic and solvation properties of electrolyte solutions with the statistical associating fluid theory for potentials of variable range, Mol. Phys. 112 (17) (2014) 2339-2364, https://doi.org/10.1080/00268976.2014.910316.

[29] M.A. Selam, I.G. Economou, M. Castier, A thermodynamic model for strong aqueous electrolytes based on the eSAFT-VR Mie equation of state, Fluid. Phase Equilib. 464 (2018) 47-63, https://doi.org/10.1016/j.fluid.2018.02.018.

[30] N. Novak, G.M. Kontogeorgis, M. Castier, I.G. Economou, Modeling of Gas Solubility in Aqueous Electrolyte Solutions with the eSAFT-VR Mie Equation of State, Ind. Eng. Chem. Res. 60 (42) (2021) 15327-15342, https://doi.org/10.1021/acs.iecr.0c02923.

[31] I.K. Nikolaidis, N. Novak, G.M. Kontogeorgis, I.G. Economou, Rigorous Phase Equilibrium Calculation Methods for Strong Electrolyte Solutions: the Isothermal Flash, Fluid. Phase Equilib. 558 (2022) 113441, https://doi.org/10.1016/j.fluid.2022.113441.

[32] N. Novak, G.M. Kontogeorgis, M. Castier, I.G. Economou, Extension of the eSAFT-VR Mie equation of state from aqueous to non-aqueous electrolyte solutions, Fluid. Phase Equilib. 565 (2023) 113618, https://doi.org/10.1016/j.fluid.2022.113618.

[33] Fufang Yang, Georgios M. Kontogeorgis, Jean-Charles de Hemptinne, Composition-dependence of relative static permittivity in ePPC-SAFT for mixed-solvent alkali halides, Fluid. Phase Equilib. 583 (2024) 114103.

[34] Fufang Yang, Abtin Raepispour Shirazi, Juan Sebastian Roa Pinto, Georgios M. Kontogeorgis, Nicolas Ferrando, Amparo Galindo, Jean-Charles de Hemptinne, Ion pairing in ePPC-SAFT for aqueous and mixed-solvent alkali halide solutions, Ind. Eng. Chem. Res. 63 (2024).

[35] C. Held, A. Prinz, V. Wallmeyer, G. Sadowski, Measuring and modeling alcohol/
salt systems, Chem. Eng. Sci. 68 (1) (2012) 328-339, https://doi.org/10.1016/j.ces.2011.09.040.

[36] M. Bülow, M. Ascani, C. Held, ePC-SAFT advanced – Part II: application to salt solubility in ionic and organic solvents and the impact of ion pairing, Fluid. Phase Equilib. 537 (2021) 112989, https://doi.org/10.1016/j.fluid.2021.112989.

[37] M. Bülow, M. Ascani, C. Held, ePC-SAFT advanced - Part I: physical meaning of including a concentration-dependent dielectric constant in the born term and in the Debye-Hückel theory, Fluid. Phase Equilib. 535 (2021) 112967, https://doi.org/10.1016/j.fluid.2021.112967.

[38] F. Yang, T.D. Ngo, J.S. Roa Pinto, G.M. Kontogeorgis, J.-C.de Hemptinne, Systematic evaluation of parameterization approaches for the ePPC-SAFT model for aqueous alkali halide solutions, Fluid. Phase Equilib. 570 (2023) 113778, https://doi.org/10.1016/j.fluid.2023.113778.

[39] S. Naseri Boroujeni, B. Maribo-Mogensen, X. Liang, G.M Kontogeorgis, Theoretical and practical investigation of ion-ion association in electrolyte solutions, J. Chem. Phys. 160 (2024) 15, https://doi.org/10.1063/5.0198308.

[40] J. Bjerrum, Adv. Chem. Ser. 62 (1966) 178.

[41] F. Yang, G.M. Kontogeorgis, J.-C.de. Hemptinne, Systematic evaluation of parameterization approaches for the ePPC-SAFT model for aqueous alkali halide solutions. II. Alkali bromides, iodides, fluorides, and lithium halides, Fluid. Phase Equilib. 573 (2023) 113853, https://doi.org/10.1016/j.fluid.2023.113853.
[42] A.S. Wexler, A step-wise ion hydration model of aqueous electrolyte solution: the 1: 1 punch, Fluid. Phase Equilib. 559 (2022) 113498, https://doi.org/10.1016/j.fluid.2022.113498.

[43] J.K. Percus, G. Yevick, Phys. Rev. 110 (1966) 251.

[44] E. Waisman, J.L. Lebowitz, Mean spherical model integral equation for charged hard spheres i. method of solution, J. Chem. Phys. 56 (6) (1972) 3086-3093, https://doi.org/10.1063/1.1677644.

[45] L. Blum, Primitive Electrolytes in the Mean Spherical Approximation, in Theoretical Chemistry, Elsevier, 1980, pp. 1-66.

[46] O. Bernard, L.B. Binding mean spherical approximation for pairing ions: an exponential approximation and thermodynamics, J. Chem. Phys. 104 (1996).

[47] L. Blum, O. Bernard, The general solution of the binding mean spherical approximation for pairing ions, J. Stat. Phys. 79 (3-4) (1995) 569-583, https://doi.org/10.1007/BF02184871.

[48] The general solution of the binding mean spherical approximation for pairing ions.

[49] M.K. Dixit, B.L. Tembe, Potentials of mean force of sodium chloride ion pair in dimethyl sulfoxide-methanol mixtures, J. Mol. Liq. 178 (2013) 78-83, https://doi.org/10.1016/j.molliq.2012.09.026.

[50] M. Kelley, A. Donley, S. Clark, A. Clark, Structure and dynamics of nacl ion pairing in solutions of water and methanol, J. Phys. Chem. B 119 (51) (2015) 15652-15661, https://doi.org/10.1021/acs.jpcb.5b07492.

[51] J.-P. Simonin, O. Bernard, L. Blum, Ionic solutions in the binding mean spherical approximation: thermodynamic properties of mixtures of associating electrolytes, J. Phys. Chem. B 103 (4) (1999) 699-704, https://doi.org/10.1021/jp9833000.

[52] A. Ruas, O. Bernard, B. Caniffi, J.-P. Simonin, P. Turq, L. Blum, P. Moisy, Uranyl (VI) nitrate salts: modeling thermodynamic properties using the binding mean spherical approximation theory and determination of "fictive" binary data, J. Phys. Chem. B 110 (7) (2006) 3435-3443, https://doi.org/10.1021/jp056733t.

[53] Pierre Turq, Lesser Blum, Olivier Bernard, Werner Kunz, Conductance in associated electrolytes using the mean spherical approximation, J. Phys. Chem. B (1995).

[54] J. Rozmus, J.-C.de Hemptinne, A. Galindo, S. Dufal, P Mougin, Modeling of strong electrolytes with ePPC-SAFT up to High Temperatures, Ind. Eng. Chem. Res. 52 (29) (2013) 9979-9994, https://doi.org/10.1021/ie303527j.

[55] P.K. JOG, S.G. Sauer, J. Blaesing, W.G. Chapman, Application of dipolar chain theory to the phase behavior of polar fluids and mixtures, Ind. Eng. Chem. Res. 40 (21) (2001) 4641-4648, https://doi.org/10.1021/ie010264%2B.

[56] K.E. Gubbins, C.H. Twu, Thermodynamics of polyatomic fluid mixtures—I theory, Chem. Eng. Sci. 33 (7) (1978) 863-878, https://doi.org/10.1016/0009-2509(78) 85176-8.

[57] J.P. Simonin, O. Bernard, L. Blum, Description of the thermodynamic properties of aqueous ionic solutions within the mean spherical approximation, Oil Gas Sci. Technol. - Revue de l'IFP 63 (3) (2008) 321-327, https://doi.org/10.2516/ogst: 2008012.

[58] J.-P. Simonin, O. Bernard, L. Blum, Real Ionic Solutions in the Mean Spherical Approximation. 3. Osmotic and Activity Coefficients for Associating Electrolytes in the Primitive Model, J. Phys. Chem. B 102 (22) (1998) 4411-4417, https://doi.org/10.1021/jp9732423.

[59] P.K. JOG, W.G. Chapman, Application of Wertheim's thermodynamic perturbation theory to dipolar hard sphere chains, Mol. Phys. 97 (3) (1999) 307-319, https://doi.org/10.1080/00268979909482832.

[60] J.P. Simonin, A. Raeispour Shirazi, O. Bernard, J.C. de Hemptinne. Simplified versions of the Mean Spherical Approximation (MSA) for the prediction of activity coefficients in electrolytes.

[61] J.S. Roa Pinto, N. Ferrando, J.-C.de Hemptinne, A. Galindo, Temperature dependence and short-range electrolytic interactions within the e-PPC-SAFT framework, Fluid. Phase Equilib. 560 (2022) 113486, https://doi.org/10.1016/j.fluid.2022.113486.

[62] F. Yang, T.D. Ngo, G.M. Kontogeorgis, Hemptinne J.-C. de, A Benchmark database for mixed-solvent electrolyte solutions: consistency analysis using E-NRTL, Ind. Eng. Chem. Res. 61 (42) (2022) 15576-15593, https://doi.org/10.1021/acs.iecr.2c00059.

[63] F. Yang, J. Qu, G.M. Kontogeorgis, Hemptinne J.-C. de, Reference Density Database for 20 Aqueous Alkali Halide Solutions, J. Phys. Chem. Ref. Data 51 (4) (2022), https://doi.org/10.1063/5.0124173.

[64] Y.-X. Zuo, W. Fürst, Prediction of vapor pressure for nonaqueous electrolyte solutions using an electrolyte equation of state, Fluid. Phase Equilib. 138 (1-2) (1997) 87-104, https://doi.org/10.1016/S0378-3812(97)00145-3.

[65] A.A. Maryott, E. S, Table of Dielectric Constants of Pure Liquids. Government Printing, NBS Circular 514 (1951). U.S.,

[66] Y.Y. Akhadov, Dielectric Properties of Binary Solutions, Pergamon, Oxford, 1981.

[67] C. Reichardt, Ion Solvation. Von Y. Marcus. John Wiley, Chichester 1985. XIII, 306 S., geb. £ 42.00. — ISBN 0-471-90756-1, Angew. Chemie 99 (8) (1987) 826, https://doi.org/10.1002/ange.19870990839.

[68] M.S. Wertheim, Fluids with highly directional attractive forces. IV. Equilibrium polymerization, J. Stat. Phys. 42 (3-4) (1986) 477-492, https://doi.org/10.1007/BF01127722.

[69] Y. Marcus, Thermodynamics of solvation of ions. Part 5.—Gibbs free energy of hydration at 298.15 K, J. Chem. Soc., Faraday Trans. 87 (18) (1991) 2995-2999, https://doi.org/10.1039/FT9918702995.

[70] M.S. Wertheim, Fluids with highly directional attractive forces.: I. Statistical thermodynamics, J. Stat. Phys. 35 (1984) 19-34.

[71] M.S. Wertheim, Fluids with highly directional attractive forces.: II. Thermodynamic perturbation theory and integral equations, J. Stat. Phys. 35 (1984) 35-47.

[72] M.S. Wertheim, Fluids with highly directional attractive forces.: III. Multiple attraction sites, J. Stat. Phys. 42 (1986) 459-476.

[73] S. Ahmed, N. Ferrando, J.-C.de Hemptinne, J.-P. Simonin, O. Bernard, O Baudouin, A New PC-SAFT Model for Pure Water, Water-Hydrocarbons, and Water-Oxygenates Systems and Subsequent Modeling of VLE, VLLE, and LLE, J. Chem. Eng. Data 12 (61) (2016) 4178-4190, https://doi.org/10.1021/acs.jced.6b00565.

[74] L. Pauling, The Nature of the Chemical Bond, Cornell University Press., 1960.

[75] S.O. Derawi, G.M. Kontogeorgis, M.L. Michelsen, E.H. Stenby, Extension of the cubic-plus-association equation of state to glycol-water cross-associating systems, Ind. Eng. Chem. Res. 42 (7) (2003) 1470-1477, https://doi.org/10.1021/ie0206103.

[76] W.R. Fawcett, Thermodynamic parameters for the solvation of monatomic ions in water, J. Phys. Chem. B 103 (50) (1999) 11181-11185, https://doi.org/10.1021/jp991802n.

[77] Yan W.D., Xu Y.J., Han S.J. Activity coefficients of sodium chloride in methanol - water mixed solvents at 298.15 K, 1994, 52, 937-946.

[78] J. Jaramillo-Gutiérrez, J. Torres-Arenas, J.-P. Simonin, O. Bernard, Description of electrolytes at high temperatures within the binding mean spherical approximation (BiMSA), Mol. Phys. (2023), https://doi.org/10.1080/00268976.2023.2287686.

[79] P.C. Ho, D.A. Palmer, R.E. Mesmer, Electrical conductivity measurements of aqueous sodium chloride solutions to 600C and 300 MPa, J. Solution. Chem. 23 (9) (1994) 997-1018, https://doi.org/10.1007/BF00974100.

[80] P. Fournier, E.H. Oelkers, R. Gout, G. Pokrovski, Experimental determination of aqueous sodium-acetate dissociation constants at temperatures from 20 to 240°C, Chem. Geol. 151 (1-4) (1998) 69-84, https://doi.org/10.1016/S0009-2541(98) 00071-0.

[81] K.D. Collins, Why continuum electrostatics theories cannot explain biological structure, polyelectrolytes or ionic strength effects in ion-protein interactions, Biophys. Chem. 167 (2012) 43-59, https://doi.org/10.1016/j.bpc.2012.04.002.

[82] M. Arrad, A.S. Wexler, A step-wise ion hydration model of aqueous electrolyte solution: the 2:2, 2:1 and 1:2 punches, Fluid Phase Equilib. 566 (2023) 113694, https://doi.org/10.1016/j.fluid.2022.113694.

[83] A.S. Wexler, Raoult Was Right After All, ACS. Omega 4 (7) (2019) 12848-12852, https://doi.org/10.1021/acsomega.9b01707.

[84] A.S. Wexler, Raoult was right after all: statistical mechanics derivation and volumetric validation, Fluid. Phase Equilib. 531 (2021) 112899, https://doi.org/10.1016/j.fluid.2020.112899.

[85] K.D. Collins, Charge density-dependent strength of hydration and biological structure, Biophys. J. 72 (1) (1997) 65-76.

[86] J.E. Combariza, N.R. Kestner, J. Jortner, Energy-structure relationships for microscopic solvation of anions in water clusters, J. Chem. Phys. 100 (4) (1994) 2851-2864, https://doi.org/10.1063/1.467231.