Importance of the Relative Static Permittivity in electrolyte SAFT-VR Mie Equations of State

The inﬂuence and importance of the relative static permittivity (RSP) in electrolyte equations of state is examined for the case of aqueous sodium chloride. Using the SAFT-VR Mie model, the Debye-Hückel (DH) or Mean-Spherical Approximation (MSA) terms, as well as the Born-solvation term, are used to formulate an electrolyte equation of state. The RSP is obtained from a variety of models, each differing in their dependencies; we consider constant, temperature, density-and composition-dependent models. For a fair comparison between different combinations of electrostatic and RSP models, all ion-related parameters are obtained a priori . A novel combining rule is proposed to obtain the unlike parameters between solvents and ions; its reliability is examined for a variety of electrolyte systems. We also compare its performance relative to parameterised electrolyte models. Both the DH and MSA terms yield similar re-sults for almost all properties and conditions. The RSP models used have the more-signiﬁcant impact. Liquid densities and solvent saturation pressures showed limited changes between RSP models whereas osmotic coeﬃcients, mean ionic activity coeﬃcients and carbon dioxide solubilities observed drastically different behaviour. Analysing the contributions of the various terms to the activities of each species in an electrolyte mixture reveals an important balance between the Born-solvation and the DH or MSA terms where the RSP models have a signiﬁcant inﬂuence over this balance, particularly when these carry a solvent-or ion-composition dependence. © 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Electrolyte systems are of growing importance in many areas of research and industry; examples include carbon capture [1] , acid gas scrubbing [2] , carriers of active pharmaceutical ingredients and formulation additives [3] , and PUREX processes [4] .To model the thermodynamic properties of such systems, one needs to account for a wide variety of interactions: short-range london-dispersion, highly directional dipolar, quadrupolar and hydrogen bonding interactions, and long-range electrostatic interactions.For each of the aforementioned interactions, various models have been developed to account for them.Short-range dispersive interactions can be modelled accurately using cubic equations of state, such as Peng-Robinson (PR) [5] or Soave-Redlich-Kwong (SRK) [6] .Whilst polar interactions and hydrogen-bonding interactions can be modelled as associative type interactions, for which the various Statistical Associating Fluid Theory (SAFT) equations of state [7][8][9][10] have proven to give excellent descriptions of, the former can be modelled as its own distinct interaction [11] .For the long-range electrostatic interactions, the two most-common approaches are the Debye-Hückel [12] (DH) and the primitive Mean-spherical approximation [13] (MSA) equations.
Although the merits of the underlying models used to characterise the dispersive and associative interactions have been the subject of various studies [37,38] , within the context of electrolyte systems, the aforementioned models typically differ in two key aspects: the model used to characterise the electrostatic interactions (DH vs MSA) and the model used to determine the relative static permittivity (RSP) of the medium.There are many studies that compared the former two approaches [22,39,40] ; one of these by Maribo-Mogensen et al. [41] performed a comparison between the DH and MSA equations, examining the behaviour of the screening length, derivatives of the Helmholtz free energy with respect to volume, temperature and composition, and the effects of the underlying parameters.It was found that the DH and MSA equations give very similar estimates of the various properties examined.However, the latter study did not examine the impact of these equations when they are part of a 'complete' model and how the real properties of electrolyte systems are affected.In addition, the comparison was done for isothermal-isochoric systems, where the impact of the electrolyte on the density is neglected.As a result, it is indeed worthwhile to re-visit this comparison, but now examining the impact on real properties and compare the predicted behaviour to experimental data.
The model providing the description of the dispersive and associative interactions is chosen to be the SAFT-VR Mie equation of state developed by Lafitte et al. [10] .Although it is likely that any of the SAFT-based models could serve the same role (allowing the conclusions of this study to be extended to other SAFT-based electrolyte models), the primary reason for the selection of this model is due to the fact that two variants have been developed which extend the model to electrolyte systems: SAFT-VRE Mie developed by Eriksen et al. [31] and eSAFT-VR Mie developed by Selam et al. [32] .The primary distinction of these models is the approach used to characterise the electrostatic interactions: the former uses the MSA equation whilst the latter uses the DH equation.Unfortunately, although a comparison between these models would be informative, there are other significant differences between them that would make a like-for-like comparison difficult.
One such difference, and one that often distinguishes electrolyte model approaches, is the choice of model used to determine the RSP.This choice is typically the source of great debate within the electrolyte modelling community, particularly with respect to what the RSP should be dependent on.Some approach use only a temperature dependence [28,32,42] , a temperature and density dependence [21,22,31,43] or a temperature, density and composition dependence [25,29] .However, no extensive comparison has been made between the difference choices of RSP models, considering its impact on the different electrostatic interaction models and the thermodynamic properties obtained.
Within this study, using the system of aqueous sodium chloride as a case study, we compare the impact of the combination of different electrostatic interaction and RSP models within a 'complete' electrolyte equation of state.To allow for a fair comparison, and to isolate the impact of the aforementioned models, we determine all model parameters a priori from either experimental data [44,45] or combining rules.In the case of the latter, we derive a novel combining rule developed to obtain the unlike interaction parameters between charged and neutral species; we provide a brief evaluation of the validity of this combining rule for a range of electrolyte systems.The effect of this parameterisation on the accuracy of the properties predicted is also considered by comparing the performance to the regressed electrolyte SAFT-VR Mie models of Eriksen et al. and Selam et al. .
In order to better understand the observed behaviour, we also analyse the magnitude of each of the contributions from the overall equation of state to the activity of each species in the electrolyte system.This also provides a better understanding of the role of the Born-solvation term within an electrolyte equation of state, which is a contribution that various electrolyte models either include [21,22,24,25,29,31,43] or ignore [26][27][28]42] .This discussion is in a similar vein to the work by Sun et al. [46] , now focusing on the SAFT-VR Mie equation, the impact of various RSP models on these contributions and each of the species present in an electrolyte system.
The subsequent sections of this article are set out as follows.Within section 2.1 , we provide an overview of the SAFT-VR Mie equation of state, the Born solvation, the Debye-Hückel and Meanspherical approximation terms; we also identify aspects in which Eriksen et al. and Selam et al. differ, beyond the electrostatic interaction model.We then give a similar overview of the different RSP models compared in section 2.2 , highlighting the different variables and parameters which they depend on.Subsequently, we illustrate how the model parameters are obtained within this study, providing a derivation for a novel combining rule between charged and neutral species ( section 2.3 ).The methods used to determine the thermodynamic properties of electrolyte systems is described in sections 2.4 and 2.5 .The results obtained from these models are given in section 3 where we examine the predicted RSP, density, isobaric heat capacity, speed of sound, osmotic coefficient, mean ionic activity coefficient, solvent saturation pressure and solubility of carbon dioxide in the aqueous sodium chloride system.Within section 4 , we analyse in greater depth the role of the various terms in the electrolyte equation of state (4.1) and the efficacy of the proposed combining rule (4.2) .Finally, in section 5 , we summarise our results and provide concluding remarks.

Theory
This section will first outline the theory behind modelling electrolytes within the SAFT-VR Mie framework.Subsequently, the various relative static permittivity models implemented within this framework will also be discussed.The means by which the model parameters for the ions are obtained will be described in section 2.3 .

Electrolyte Equation of State
As with any system in equilibrium, one can often rely on equations of state to obtain a wide variety of thermodynamic properties.Typically, these equations can be split between an ideal and residual contribution where, if expressed in terms of the Helmholtz free energy, A : where N, k B and T denote the number of particles in the system, the Boltzmann constant and temperature, respectively.The subscripts ideal and res.denote the ideal and residual properties, respectively.In the case of electrolyte systems, the residual contributions can account for a variety of interactions: short-range dispersive interactions, highly-directional associative interactions ( e.g.hydrogen bonding), and long-range electrostatic interactions.The first two of these interactions has been very effectively captured by the many variants of the SAFT equation of state [9,10] .Whilst many of these variants have been applied to electrolyte systems [25][26][27][28] , the one that will be of interest in this study is the SAFT-VR Mie equation developed by Lafitte et al. [10] .
Another reason for the selection of the SAFT-VR Mie equation in this study relates to how one can account for electrostatic interactions of electrolyte systems within this framework.To date,  [31] and eSAFT-VR Mie developed by Selam et al. [32] .Whilst these models differ in more than a few elements, the primary distinction between the two is the choice of model to account for electrostatic interactions between ions.Although both models implement the same Born solvation term, the former uses the MSA model [13] and the latter opts for the DH model [12] .The debate as to which, between the MSA and DH models, is best suited to model electrolyte systems has been lengthy; Maribo-Mogensen et al. [41] have shown that, up to a scaling of the ion parameters, the two approaches give similar results.
As a result, along with the Born solvation term, both models are implemented and compared within this study.

Ideal Contribution
Often taken for granted, the ideal contribution to any equation of state is of vital importance when evaluating thermodynamic properties, even in the liquid phase [47] .It can generally be expressed as a sum of the translation, rotational and vibrational contributions: where x i , ρ i , i , n rot ,i and θ rot ,i denote the molar composition, number density, thermal de Broglie wavelength, number of rotations and rotational temperature of species i .n vib ,i , g i, v and θ vib ,i, v denote the number of vibrational modes on species i and, the vibrational degeneracy and temperature of mode v on species i , respectively.n species denotes the number of species present in the system.As implemented by Walker and Haslam [47] , the rotational temperature for all species is set to 1 K.For monoatomic ions, the number of rotations and vibrations is set to zero; as the vibrational temperatures for polyatomic ions are not currently available, properties involving this contribution will be avoided.For water and carbon dioxide, the parameters are given in Table 1 .

SAFT-VR Mie Contribution
As introduced by Lafitte et al. [10] , species are modelled as chains of spherical segments of length m with short-range dispersive interactions characterised by an underlying Mie potential: where r is the inter-segment distance.When i = j, , σ , λ r and λ a are the potential well depth, segment diameter and, repulsive and attractive exponents of the segments in a given species, respectively.For i = j, the parameters characterise the unlike interactions between the segments of two species.The variable C i j is given by: which ensures that the minimum of the potential is alwaysi j .In addition, highly directional interactions, such as hydrogen bonding, are represented using association sites present on the segments whose interactions are modelled using a square-well potential: where r ab is the center-center distance between association sites a and b.HB ab and r c ab are the potential well depth and width of the square-well potential, respectively.
Using these parameters, for a given composition ( N ), temperature ( T ) and volume ( V ), the residual Helmholtz free energy due to the interactions described above can be obtained by summing over three contributions: where the subscripts indicate the particular contribution to the residual Helmholtz free energy: mono.denotes the contribution from the presence of the spherical segments and dispersive interactions between segments, chain denotes the contribution from the formation of a chain of the Mie segments and, assoc.denotes the contributions from the interactions between association sites present on the segments.
In passing, we note that the Barker-Henderson theory [4 8,4 9] plays an integral role in the development of SAFT-VR Mie equation.As such, the Barker-Henderson hard-sphere diameter is a key variable within this equation; it is evaluated using the following equation: where β = 1 / (k B T ) .As it is not possible to obtain an analytical solution to this integral, numerical methods are required.Lafitte et al. [10] suggest that the use of 10-point Gauss-Legendre quadrature should give a sufficiently accurate approximation of the integral [50] .However, more recent implementations (particularly in SAFTγ Mie [51] ) use 5-point Gauss-Laguerre quadrature instead [52] .
As shown in Fig. 1 a, whilst the latter approach is more accurate in a lower reduced temperature range ( T * = T / ), for ions can be of the order of magnitude of 10 0 .Thus, at room temperature, T * ∼ 100 where Gauss-Legendre quadrature provides a better approximation (see Fig. 1 b).As such, careful considerations must be made when deciding which quadrature to use (we point out that both Eriksen et al. [31] and Selam et al. [32] have used 10-point Gauss-Legendre quadrature).
We now examine the association term as it is of particular importance for electrolytes, especially in the context of the RSP.It is given by: where X a,i is the fraction of association sites a on species i that are not bonded to any other sites, n sites ,i is the number of site types on species i and n a,i is the number of sites a on species i .It is evaluated by solving a system of mass-action equations: b=1 n b, j X b, j ab,i j (9) where ab,i j characterises the association strength between sites a and b on species i and j, respectively.Within the SAFT-VR Mie framework, there are three ways one can determine ab,i j , the first of which, as given by Lafitte et al. [10] uses the following equation:  where F ab,i j is the Mayer function of the bonding interaction, g HS d,i j (d) is the pair-distribution function for a hard-sphere fluid of diameter d (as denoted by the subscript) and κ ab,i j is the so-called dimensionless bonding volume.However, Dufal et al. [53] proposed two alternatives, both of whom can be expressed as: ab,i j = F ab,i j I i j K ab,i j , (11) where I i j is referred to as the association kernel and K ab,i j is also bonding volume which now carries units of volume.Two such kernels were developed: one referred to as the Lennard-Jones (LJ) kernel (where I i j (σ i j , i j ) ) and the other referred to as the Mie kernel (where I i j (σ i j , i j , λ r ,i j ) ); along with these kernels, two water models were parameterised.The one obtained from the LJ kernel became the 'default' in most SAFT-VR Mie (and SAFT-γ Mie) models and, as such, binary interaction parameters between water and other species are readily available in literature [54,55] ; this is indeed the model used by Eriksen et al. .However, the model obtained from the Mie kernel is the one used by Selam et al. .The SAFT-VR Mie parameters for non-ionic species used within this study are given in Table 2 .How the parameters for ionic species are obtained is be outlined in section 2.3 .

Born Solvation Term
The Born solvation term accounts for the effects of solvation when an ion is introduced to a solvent / medium of relative static permittivity r .The Gibbs free energy change for taking an ion of diameter σ Born i from a vacuum to the medium (equivalent to the Helmholtz free energy when the relative static permittivity is constant) is: where e and 0 are the elementary charge of an electron and permittivity of the vacuum, respectively.Z i is the charge of i and σ Born i is the effective radius of ion i .The physical interpretation of this radius is the topic of much debate; in some frameworks, it is equivalent to σ in the SAFT model [42,43] , whereas others treat it as a separate radius related to the cavity occupied by the ion in the solvent.The latter interpretation is the one adopted by Eriksen et al. and Selam et al.where both obtain this radius from experimental values [44] .Models developed in this work will also use these values for σ Born i .

Debye-Hückel Term
One of the first models to characterise the electrostatic interactions, the DH equation [12,57] can be derived by assuming that the distribution of ions follows a Boltzmann distribution and solving the linearised Poisson equation to give us the electrical potential.From this, the Helmholtz free energy is given as: where κ is the inverse screening length given by: and the function χ is expressed as: (1 where: It should be noted that in

Mean-Spherical Approximation Term
Along with the DH model, the primitive mean-spherical approximation model developed by Blum [13] is another approach often used to characterise the electrostatic interactions within a solvent.The primary difference between the two approaches is whether the ions are treated as hard-spheres, leading to an excluded volume inaccessible to other ions.The Helmholtz free energy within this framework is given by: where the excess internal energy contribution, U MSA , is given by: The screening length, , can be determined from: where the effective charge, Q i , is given by: and the auxilliary functions and P n are given by: The packing fraction, , has been defined one of two ways.Within the SAFT-VRE Mie model (and others [42,43] ), it is given by: whereas, Maribo-Mogensen et al. [41] express it as: The only distinction between these two equations is that, in the latter, the solvent is included within the sum giving it a dependence on the solvent concentration.It is not readily apparent in Blum [13] 's original work as to which one is correct, however, examining the impact of this dependence on the chemical potential of the solvent in Maribo-Mogensen et al. [41] 's work, this dependence is orders of magnitude smaller than it is for the ions.
As such, and to remain comparable to Eriksen et al. 's model, Eq. ( 23) will be used instead.Note that the screening length is dependent on itself and must therefore be solved iteratively where the Debye screening length (Eq.( 14) ) serves as a good initial guess.
We take this opportunity to point out the existence of the nonprimitive MSA (NP-MSA) model, also developed by Blum and coworkers [33,34] .This approach, in contrast to the primitive approach, explicitly accounts for the structure of the solvent, including ion-ion, ion-dipole and dipole-dipole interactions.One interesting aspect of this approach is, as observed by Simonin [58] , due to the correct description of the ion-solvent interactions, the 'Born' energy is already accounted for.Although the SAFT-VR Mie framework has not yet been extended to include the NP-MSA model, the SAFT-VR SW framework has in the form of the SAFT-VR+DE equation of state [35,[59][60][61] .This more-accurate representation of the solvent has allowed the SAFT-VR+DE equation to obtain results

Table 3
Parameters for Eq. ( 25) corresponding to water.comparable or superior to those obtained from both SAFT-VR Mie electrolyte models, although the former regressed salt-specific parameters, rather than ion-specific parameter in the case of the latter models.
Another interesting aspect of this model is that the relative static permittivity used is found to be dependent on the composition of both the solvent and ions [33,62] .As this hints towards the need for an ion-composition-dependent RSP model, it will be interesting to observe how including this dependence will impact the performance of the electrolyte models.

Relative Static Permittivity models
Another facet we intend to explore which, in the context of electrolyte SAFT-VR Mie equations of state, has not been explored in great depth in previous studies is the role of the relative static permittivity within this model.How one should model the RSP within the context of either MSA or DH theory is a topic of debate.Naturally, the most straight-forward approach is to keep it constant; this shall indeed be considered in this work with a value of r = 78 .4 (RSP of water at 298 K).However, works such as Selam et al. have included a temperature dependence to the RSP.Others, such as Eriksen et al. [31] , have also included a density and electrolyte-free concentration dependence of the solvent in their models of the RSP.Finally, one could also consider the dependence of the RSP on the concentration of the ions although such approaches are difficult to discuss in primitive models such as MSA and DH theory.
All of these approaches to modelling the RSP shall be considered within this work.The models were selected based on varying levels of complexities i.e. they account for the temperature, density and composition (solvent and electrolyte) dependence of the RSP.

Zuo and Fürst's correlation
The correlation used within eSAFT-VR Mie model, Zuo and Fürst [63] 's accounts for the temperature dependence of the RSP of water, parameterised using five θ i values given in Table 3 with the equation given below: The parameters are fitted within a temperature range of 288.15 K-403.15K. We note that the above correlation does equal zero at approximately 600 K which will lead to singularities when evaluating properties in this temperature range due to the presence of −1 r in various equations.

Valiskó and Boda's correlation
This correlation developed by Valiskó and Boda [64] and only depends on the concentration of the electrolyte; for example, in the case of aqueous sodium chloride: where c is the concentration of the sodium chloride (NaCl) in mol dm −3 .This correlation was fitted in range of 0 mol kg −1 to approximately 5 mol kg −1 .However, it is limited to a temperature of 298K and pressure of 1 atm.It is worth pointing out that this correlation also has a dependence on the density and, as a result, the accuracy of this correlation also depends on the accuracy of the density estimated from our equation of state.Nevertheless, many such correlations exist for other salts and have been implemented within both DH and MSA theories [65][66][67] .

Schreckenberg et al. 's correlation
Originally developed by Schreckenberg et al. [43] for the SAFT-VRE SW equation [43] , this model accounts for the temperature and density dependence of the RSP, along with the dependence on the solvent concentration: where ρ solvent denotes the number density of the solvent and d is obtained from: where x 0 ,i denotes the electrolyte-free solvent composition and d i j is a solvent-dependent parameter where, for i = j: where d V,i and d T,i are also solvent-specific parameters; for species of interest in this study, these have been tabulated in Table 4 .For i = j, d i j is simply obtained from the arithemic mean: The parameters for water were obtained by fitting to experimental data in a temperature range of 273 K to 423 K and densities up to 63 mol dm −3 .

Michelsen and Mollerup's correlation
Michelsen and Mollerup [68] suggested that the effect of the presence of ions within a solvent could be modelled through a correction to the pure RSP of the solvent in the following way: (31) where E(V, N ion ) is the ion correction factor.The RSP of water, in this model, is obtained from: where T 0 = 273 .15 K is the reference temperature, 0 r ,w = 87 .82 is the RSP of water at T 0 , μ w, 0 = 1 .855 D is the dipole moment of water, N A is Avogadro's constant and β 1 = 3 .1306 is a fitting constant.The water density ρ w is obtained from the volume obtained from the equation of state.The ion correction factor can be expressed as: where c k is the molar concentration of ion k in mol L −1 , β 2 = 0 .160 L mol −1 , β 3 = 0 .010 L mol −1 and α k is an ion-specific constant where α Na + = 0 .1062 L mol −1 and α Cl − = 0 .1172 L mol −1 .This correlation is compatible with more ions than the ones of interest in this study.It also provides a temperature and, albeit minor, volume dependence of the RSP.
This correlation has been used by Maribo-Mogensen et al. [41] in a comparison of the DH and MSA theories in modelling electrolyte properties.

Maribo-Mogensen et al. 's models
Taking a more physically-based approach, Maribo-Mogensen et al. [69] start from Hasted [70] 's expression for the RSP: where g i is the Kirkwood g-factor introduced by Fröhlich and Maradudin [71] and ∞ is the RSP at infinite frequency which can be obtained from the Clausius-Mossoti equation [72,73] : where α 0 ,i is the molecular polarizability of species i (given in Table 5 ).The Kirkwood g-factor is obtained from the following equation: where z i j is the coordination number of species j around central molecule i , γ i j is the angle between the dipole moments of molecules i and j in the hydrogen bonding network ( γ ii = 69 .4 deg for water), θ i j is the hydrogen bond angle between molecule j in the shell around molecule i and the second shell neighbour ( θ ii = 108 .9 deg for water).Note that these last two parameters are typically fitted to experimental data.P i j is the probability molecule i is associated to molecule j, determined by the following set of equations P ab,i j = ρx j X a,i X b, j ab,i j , (38) and P i is the probability of molecule i being associated with any other molecule, given by: In the above equations, X a,i and ab,i j are both the same variables one obtains from the association term in SAFT-VR Mie (Eq.( 9) and ( 10) ).We note here that the parameters γ i j and θ i j are fitted to the experimental values of the RSP where the association fractions are obtained from the cubic-plus association (CPA) equation of state [23] .As such, this model is effectively predicting when implemented with the SAFT-VR Mie equation of state.However, this model ignores the short-range ion-dipole interactions.If one wants to examine the influence of the presence of ions on the RSP, one could still use Eq. ( 34) if ion-solvent and ionion association were explicitely accounted for within the association term of the SAFT-VR Mie equation; in most models, this is not the case (see Gil-Villegas et al. [74] for a detailed discussion on the impact of including these interactions in an electrolyte model).As such, Maribo-Mogensen et al. [75] developed a modification to Eq. ( 34) that would account for this association: where i is the fraction of component i that is not bonded to an ion.It can be obtained by solving the following mass-action relation (where i indexes the solvent and j and k index the ions): Table 5 Various parameters needed to model the species of interest.Polarizabilities [45,77,78] , ionisation energies [77,79,80] , segment diameters [44,54,56] , Born diameters [81] ( * denotes those where the Born diameter was obtained by 1 . 07σ ), dipoles [82] and quadrupoles [83,84] are all obtained from the literature.
α 0 /( Å 3 ) where m is the molality of the salt in mol kg −1 of solvent, n i is the number of moles of species i , ν j is the stoichiometric factor of ion j in the salt and N sites , j is the number of association sites on ion j.The association strength, i j is now determined from the following equation (when the solvent in question is water): where N i j is the apparent coordination number of i around ion j.
As described by Maribo-Mogensen et al. [75] , in the case of N w j , if j is a monovalent cation, it is taken to be six, if j is a bivalent cation, it is taken to be eight, if j is an anion, it is taken to be five.
From N w j , we can obtain N sites , j : Whilst N i j between ions was not discussed explicitly within the original study, it will be taken to be the lattice coordination number of the salt.Nevertheless, this model consists of mesoscopic electrostatic equations which neglect the complex nature of the fluctuations of the dipolar moments in ion-dipole interaction length scale, particularly for small or multivalent ions [76] .

Model parameters
Whilst the primary focus of this study is aqueous sodium chloride, Table 5 provides a summary of all the parameters required to model a range of systems using the strategy developed in this section.

Ion-ion interaction parameters
For simplicity and consistency with other approaches [31,32] , all ions will be modelled as Lennard-Jones spheres ( i.e. m = 1 , λ r = 12 and λ a = 6 ).The segment size, σ , has been obtained a number of different ways.Some models [32,43] treat σ as an adjustable parameter fitted to experimental data whilst others [31]  σ used for ions will be the same.
To obtain the binary interaction parameters ( i = j) between ions, we shall use the standard combining rules for the SAFT-VR Mie equation [10,53] .For σ i j , a simple arithmetic mean is used: For the potential exponents, the combining rule is obtained by applying a geometric mean to the van der Waals attractive energy, as derived by Lafitte et al. [10] : In both Eriksen et al. and Selam et al. , the potential well-depth is obtained from a combining rule.Both only consider the London dispersion interactions between ions given by the following potential: where I i is the ionization energy of species i (given in Table 2 ).We can then obtain the integrated form of this potential from: giving: Where Eriksen et al. and Selam et al. differ is which potential one should integrate and equate to Eq. ( 50) .In SAFT-VRE Mie [31] , the full Mie potential is used, given in Eq. ( 3) , where, when integrated, gives: On the other hand, only the attractive part of the Mie potential is considered in eSAFT-VR Mie [32] : When integrated, one gets: One can then equate either 51 or 53 to Eq. ( 50) to obtain the following combining rules: These equations are essentially the Hudson and McCoubrey combining rules applied to a Mie potential [85] .Apart from when using Selam et al. 's model, Eq. ( 54) will be used to obtain the ionion like and unlike interaction parameters.However, as Selam et al. point out, it is likely that either of these combining rules will lead to similar predictions, particularly when considering that the ionion interactions are not as significant as the solvent-solvent and solvent-ion interactions at low molalities.

Ion-solvent interaction parameters
Within this section the index i shall denote ions whilst j shall denote solvents.In order to make a 'fair' comparison between the DH and MSA models (coupled with the various RSP models), it is best to use the same parameters between them such that only the models themselves are compared.Although the ion-ion interaction parameters, can be obtained a priori from the previous section, this is much more-difficult for the ion-solvent interaction parameters.
For σ i j and λ k,i j , the combining rules given by Eq. ( 46) and (47) , respectively, shall be used for the ion-solvent interaction parameters.
However, the unlike potential well depth ( i j ) is almost never obtained a priori ; indeed, it is more commonly treated as an adjustable parameter [28,31,32,43] to be fitted against experimental data.One can see why that may be the case when comparing the estimates from Eq. ( 54) to the regressed values obtained from the works of Eriksen et al. and Selam et al. in Table 6 .
It is clear that the combining rule is significantly underestimating the value of i j between ions and water when compared to regressed values.It is likely that additional interactions between the solvent and ion are present, aside from the London dispersion interactions captured in the Hudson-McCoubrey combining rule.Haslam et al. [86] have developed a combining rule that accounts for the London dispersion, dipole and quadrupole interactions present between two species modelled by the following potential: where Q 0 , j is the quadrupole moment of species i .For carbon dioxide and water, values are available in literature [82,83,87,88] , however, these can vary with temperature.Haslam et al. [86] showed that this is indeed important.However, in this work, as most of the conditions examined will be in the liquid phase at room temperature and pressure, the value of μ 0 ,i corresponding to these conditions will be used.
Integrating using Eq. ( 49) and equating the result to ψ Mie i j results in the following combining rule: Note that the number of segments m i are now included within the combining rule, in spite of not being present in Eq. ( 56) .This is because Haslam et al. [86] showed that the binary interaction potential should indeed incorporate the number of segments.When applying Eq. ( 59) , one obtains the parameters given in Table 6 .Whilst the value of i j has indeed increased, it still less than the values obtained by Eriksen et al. and Selam et al. .In addition, like the Hudson-McCoubrey combining rule, the trend between the molar mass of the cation and i j is inverted compared to the regressed values for the cations.It is very likely that an important interaction between solvents and ions has been omitted.
Most likely, this would be the ion-dipole interactions which can be captured by the following potential: θ can be eliminated from this equation by taking a rotational average: Integrating using Eq. ( 49) : This term can then be added to Eq. ( 59) to give: 3 2  54), ( 59) and ( 63) (for different values of r ).Parameters used in the combining rules are given in Table 5 and, where a temperature dependence is present, T is set to 298 K.

Ion
Eriksen et al. [  Note that r has now been included in Eq. ( 63) ; this is primarily because Haslam et al. [86] highlighted that, particularly in the liquid phase, the RSP plays an important role in obtaining accurate binary interaction parameters.This is not done in Eq. ( 59) primarily because, as shown in Table 6 , i j is already underestimated when r = 1 (further increases will only worsen this estimate).Applying Eq. ( 63) with r = 1 , the estimated i j is significantly overestimated compared to regressed values, as shown in Table 6 .This indicates that the ion-dipole term is dominant, where the remaining contributions can be neglected.Clearly the RSP must have a value other than one; the next best estimate would be the value of water at 298 K (78.4).Doing so underestimates the value of i j compared to the regressed values.Without doing any regressions, it is difficult to determine what the 'correct' value of r should be and clearly an example as to why regressing this parameter is typically required.
Without introducing a composition dependence from the models introduced in section 2.2 , a potential approximation may be a value of r = 32 , which has been used as an approximate for the RSP in salt-water systems [89,90] .Applying this to Eq. ( 63) , one obtains the values shown in Table 6 and Fig. 2 .The magnitude of It is interesting to consider, if the NP-MSA approach had been used, which accounts for the ion-dipole interactions, perhaps the Hudson-McCoubrey combining rules would be sufficient to obtain the ion-water interaction parameter as the SAFT-VR Mie terms would no longer need to account for these interactions.Unfortunately, incorporating this approach within the SAFT-VR Mie framework is beyond the scope of this study.

Thermodynamic properties of electrolytes
Electrolyte systems are typically defined using a salt molality, m , which, in this section, carries units of moles of salt per kilograms of solvent.If the salt-free composition of the solvent, x j, 0 is given, one can obtain the composition of the salt in the electrolytesolvent mixture from: where M r ,i denotes the molar mass of species i and v i is the stoichiometric coefficient of ion i in the salt.Similarly, for the ions formed from the electrolyte assuming complete dissociation in the solvent: With the conditions of the system specified ( N , T , p), one can then utilise the equation of state to obtain a variety of properties provided they solve for the corresponding volume at a given pressure p. Bulk properties such as the density, heat capacity and speed of sound are all of interest for electrolyte systems, with the latter two seldom being used to evaluate the validity of an electrolyte model.Properties that are of particular interest for electrolyte systems are the mean ionic activity coefficient ( γ ±,m , MIAC) and osmotic coefficient ( ).These can both be related to fugacity coefficients obtained from the chemical potential of a species i : where ϕ i denotes the fugacity coefficient of species i , μ i, res denotes the residual chemical potential of species i and Z denotes the compressibility factor of the system ( Z = pV Nk B T ).The activity coefficient of species i is then given by: where γ i denotes the activity coefficient of species i and the superscript denotes a property relating to a pure system of species i .However, in the case of ions, as there is not a physically meaningful pure system one can refer to, the rational asymmetric scale is used where: now the superscript * denotes a property relating to the system at infinite dilution of ion i and γ i, * now denotes the activity coefficient in the rational asymmetric scale.As x i → 0 , γ i, * → 1 .We can then convert this to a molal-based scale: where x j is the molar composition of the solvent.Taking a geometric mean (weighted by the stoichiometric coefficients of the ions in the salt) then gives the MIAC [91] : where the subscripts + and -refer to the property relating to the cation and anion in the salt, respectively.Using γ i , we can also obtain the osmotic coefficient of a solvent j: We note that the MIAC and osmotic coefficient are not independent as both must satisfy the restricted Gibbs-Duhem relation through the following relationship [92] for a single salt solution:

Phase Equilibria
For a system to be considered in thermodynamic equilibrium, it must satisfy the conditions of thermal, mechanical and chemical equilibria.For two phases in equilibrium ( α and β) these can be expressed as: However, in electrolyte systems, the conditions of chemical equilibrium are more complicated.For neutral species, chemical equilibrium can still be defined as the equality of their chemical potential in each phase: However, when the electrolyte solute is fully ionised, a different constraint must be used for the charged species across two neutral phases, as derived by Großmann and Maurer [93] : This system of equations is the one used by As we will be using other RSP models which do not account for the density of the phase, we will follow Selam et al. 's approach in constraining ions to only be present within the liquid phase.We do note that Schreckenberg et al. [43] showed that, even when implementing Eq. ( 76) , the ions would be present in the gas phase with compositions of the order of magnitude of 10 −80 , showing that for a single mole of an electrolyte system, we do not anticipate even a single ion atom to be present in the gas phase.

Results
As we aim to compare a great many combinations of RSP and electrostatic interaction models, as well as the existing SAFT-VR Mie electrolyte models [31,32] , we shall focus most of our attention on a single electrolyte system: aqueous sodium chloride.This system is selected primarily for the wealth of experimental data available which can be used to analyse the performance of each of these models.The properties we aim to cover include the RSP, density, second derivative properties, osmotic coefficients, mean ionic activity coefficients, solvent saturation pressures and solubility of carbon dioxide in the presence of an electrolyte.This should serve as an extensive comparison of all the models developed in section 2 .Throughout this section, we define the percentage average absolute deviations ( % AAD) of a given property, X, as: where i is a given data point and n data is the number of data points.The superscripts exp.and calc.denote experimental or calculated properties.For simplicity, we abbreviate each of the relative static permittivity models in the figure legends as: • (1)  r : Constant relative static permittivity (78.4) • (2)   r (T ) : Zuo and Fürst's correlation • (3)   r (ρ, x ) : Valiskó and Boda's correlation • (5)   r (T , ρ, x ) : Michelsen and Mollerup's correlation

Relative Static Permittivity
Although the permittivity is a property used within the Born solvation, DH and MSA terms, in some of the models discussed in section 2.2 , its value does depend on the liquid density predicted by the SAFT-VR Mie equation (coupled with the electrolyte-related  [72,94,95] of the relative static permittivity of an aqueous solution of sodium chloride predicted using different electrolyte SAFT-VR Mie models and RSP models.Temperature range of 273.15K to 623.15 K, pressure range of 0.1 MPa to 50 MPa and molality range of 0.0 mol kg −1 to 5.5 mol kg −1 .terms) as well as the parameters used within these equations (in the case of the Maribo-Mogensen models, it also depends on the association fraction of the solvents).As such, it would be important to consider the accuracy of the predicted permittivity between these different models.With experimental data obtained for an extensive range of conditions [72,94,95] , the % AAD of the predicted r from experimental data using different models are given in Table 7 .Note that the %AAD for the constant permittivity has also been included in spite of its obviously poor accuracy; this serves more as a 'baseline' to compare the model models.In addition, in Fig. 3 , we show the temperature ( Fig. 3 a), pressure ( Fig. 3 b) and compositional ( Fig. 3 c) dependence of the various RSP models.As it would be impractical to plot a total of 23 different models in a single figure, we examine only those which provide useful insight.

Electrolyte
At first glance of Table 7 , it is apparent that the accuracy of the RSP model is mostly independent of the electrostatic interaction model used, with the %AAD being consistent between each row.It is also the reason that we only consider the DH model in Fig. 3 .
For Zuo and Fürst [63] 's correlation (model 2), given that it does not depend on any of the parameters or properties obtained from the electrolyte model, it is not surprising that the %AADs are entirely independent of the electrostatic interaction model used.We do point out, however, that Zuo and Fürst [63] 's correlation does eventually predict a RSP with a value of zero at approximately 600 K (as shown in Fig. 3 a).Although water would only be in the liquid phase at this temperature if there is a very high pressure, this will lead to a divergence originating in the electrolyte-related terms of these models (as these are proportional to −1 r ).It is important to re-iterate that the parameters of Zuo and Fürst [63] 's correlation were not fitted to these conditions and are therefore extrapolated values.
For models such as Valiskó and Boda [64] , Schreckenberg et al. [43] and Michelsen and Mollerup [68] (models 3, 4 and 5, respectively) that have a dependence on the density, there is at most a 0.05% difference in %AAD between the different electrostatic term.The case with the largest difference is when Schreckenberg et al. [43] ' correlation (model 4) is coupled with either Eriksen et al. 's model or the models obtained through the use of Eq. ( 63) , with the former having the slightly larger %AAD.This is unexpected considering the model parameters used by Eriksen et al. have been regressed to experimental data.However, this difference is very slight.
Aside from the constant permittivity, the largest %AAD is achieved by Valiskó and Boda [64] 's correlation (model 3), which is perhaps due to the fact that it was developed to only account for the electrolyte concentration dependence of the RSP at 298 K whilst the data set used to evaluate the %AAD included mainly pure-solvent data covering an extensive range of temperature and pressure.This is partly reflected in Fig. 3 c where only the models which account for the ionic composition dependence of the RSP accurately capture experimental data; in spite of being more accurate in capturing this trend, Valiskó and Boda [64] 's correlation still has a large %AAD.We do note, however, that in spite of being designed to capture this trend, Valiskó and Boda [64] 's correlation does not capture the experimental data as well as other models in Fig. 3 c; this could be as a result of the accuracy of the volume determined by the electrolyte model, particularly when considering the work by Maribo-Mogensen et al. [25] who showed that, to obtain an accurate estimate for the density, one may require a Peneloux volume correction [97] (the subsequent discussion will examine the accuracy of the density obtained from these models).
In overall terms, the most accurate estimates of the RSP came from the Maribo-Mogensen et al. [69,75] models (6a and 6b), with the electrolyte version achieving the lowest %AAD overall.Both the electrolyte and solvent versions are identical when considering a pure solvent.As such, the improvement in %AAD between the two in Table 7 is solely due to the improved prediction of the ionic composition dependence shown in Fig. 3 c.Even when examining the temperature and pressure dependence in figures 3 a and 3 b, respectively, except when approaching room temperature and atmospheric pressure (where Schreckenberg et al. [43] 's correlation performs better) the predictions of the Maribo-Mogensen et al. [69,75] models are in excellent agreement with experimental data.
Considering their role in the electrolyte-related terms, one might be tempted to presume that the accuracy of the RSP model would be reflected in the predictions of other properties.This will be examined in subsequent sections.[96] of the liquid density of an aqueous solution of sodium chloride predicted using different electrolyte SAFT-VR Mie models and RSP models.Temperature range of 273.15K to 623.15 K, pressure range of 0.1 MPa to 10 MPa and molality range of 0.0 mol kg −1 to 5.5 mol kg −1 .

Density
The liquid density of an electrolyte system is typically a property of great interest and, as it was done by both Eriksen et al. and Selam et al. , is often used to regress the model parameters.%AAD from experimental data [96] of the different models are given in Table 8 .In Fig. 4 , the predicted liquid density are shown using different RSP models and electrolyte models, as well as comparing these to Eriksen et al. 's and Selam et al. 's models.
When examining Table 8 , between the models based on the framework discussed in section 2.3 , it is clear that the impact of the RSP model on the accuracy of the predicted densities is more significant than that of the electrostatic interaction models, with most %AADs being consistent between the DH and MSA-based models for the same RSP model, as it is in Table 7 .This is further reflected in Fig. 4 where the predictions for different electrolyte models are almost indistinguishable, particularly in the case for Zuo and Fürst [63] 's correlation (model 2).Interestingly, as we introduce the solvent-composition ( (4)  r ) and electrolyte-composition  8 ; this is reflected in Fig. 4 .This may indicate that, in spite of the improved results with a more accurate permittivity model, the σ ii parameter has a more-significant impact on the liquid density of the mixture.Indeed, other works [26,43] often treat σ ii as an adjustable parameters and obtain more-accurate liquid densities.Further comparisons could be carried out for different experimental values of σ ii from Shannon [44] 's work in a similar way to the work by Valiskó and Boda [64] for the DH term.This does indicate that if one were to regress the model parameters, not only should they include liquid densities within the data set used, but they may need to consider treating σ ii as an adjustable parameter rather than determining it a priori .

Second derivative properties
Seldom considered when evaluating an electrolyte, second derivative properties are still of interest for such systems.When these are examined, typically only the isobaric heat capacity is considered [25,100] .However, in this work, we will also consider the impact of electrolytes on the speed of sound.The %AAD from experimental isobaric heat capacities [98] and speed of sound [99] are given in Table 9 .In addition, Fig. 5 visually illustrates the impact of different electrolyte models on the isobaric heat capacity ( figures 5 a and 5 b) and speed of sound ( figures 5 c and 5 d).
The most noteworthy of these results are the %AADs for the heat capacities obtained using Zuo and Fürst [63] 's correlation (model 2) coupled with any of the electrolyte models where values are much larger than 100%.The reason for this becomes apparent when examining Fig. 5 a where a divergence towards −∞ is observed; remembering Fig. 3 a, it is clear that this divergence arises primarily as a result of the zero in the value of r .This is a rather extreme example of the importance of a reliable RSP model when extrapolating.This does indeed imply that Selam et al. 's model will be limited to certain temperatures and pressures.
Another interesting divergence is seen when using the constant permittivity model in Fig. 5 a where, as we are above the critical pressure of water, we expect to observe a divergence or maximum in the isobaric heat capacity, but this would typically occur above the critical temperature as well.Here, we see the divergence occurring at a lower temperature as we increase the molality of the electrolyte.It is expected that, as more electrolyte is introduced, the critical temperature of the mixture decreases [101,102] ; this downwards shift in the maximum is likely a reflection of this.It does raise the question as to the impact of electrolytes of supercritical liquid-gas boundaries such as the Widom-Fisher and Widom lines [103,104] , although such studies are beyond the scope of this work.
Examining the %AADs in Table 9 , it is interesting to see that the use of either the DH or MSA model had very little impact on the accuracy of the overall model, with the RSP model having the more noteworthy impact.Even in the case of Eriksen et al. 's model with regressed parameters, the %AAD is comparable to the models developed in section 2.3 .
In the case of the speed of sound, neither the RSP model or the electrostatic interaction model had a significant impact on the %AAD.From figures 5 c and 5 d, it is clear why this is the case as the underlying water model gives poor predictions of the speed of sound even when no electrolyte is present [56] (this is also the case for the isobaric heat capacity; see Fig. 5 b).However, even when one only considers the dependence of the speed of sound on the molality of the electrolyte, the trend predicted by the various %AAD from experimental data [98] of the isobaric heat capacity and speed of sound [99] of an aqueous solution of sodium chloride predicted using different electrolyte SAFT-VR Mie models and RSP models.For the isobaric heat capacity: temperature range of 323.15K to 573.15 K, pressure range of 20 MPa to 98 MPa and molality range of 0.4 mol kg −1 to 6.0 mol kg −1 .For the speed of sound: temperature range of 273.15K to 353.15K, pressure of 0.1 MPa and molality range of 0.0 mol kg −1 to 1.0 mol kg  This does indicate that, if one desires accurate second-derivatives properties, obtaining an accurate underlying solvent model is a key first step; Graham [106] 's water model might be more-well suited as it does result in more-accurate second-derivative properties [47] .However, there are discrepancies within the electrolyte models used that should also be addressed, especially for the speed of sound.

Osmotic and Mean Ionic Activity Coefficient
When evaluating the performance of an electrolyte model, two properties that are most often used as a benchmark are the mean ionic activity coefficient (MIAC) and osmotic coefficient.Both of these depend on the chemical potential of the cation, anion or solvent.Both of these properties are related through Eq. ( 72) implying that if a model performs well with one property, it should also do well in the other.It is for this reason we examine these properties simultaneously: %AAD from experimental data of the various models are provided in Table 10 , with figures 6, 7 and 8 serving as visual aids to analyse these results.
When examining Table 10 , in comparison to the properties in the previous sections, the difference in values between the models developed in section 2.3 and those of Eriksen et al. and Selam et al. is more-significant, especially for the MIAC.This perhaps indicates that, in order to obtain accurate estimates for the osmotic  coefficient and MIAC, regression is required, particularly with respect to the i −water parameter when comparing the model developed by Eriksen et al. and those developed in this work using Schreckenberg et al. [43] 's correlation (model 4), given all other parameters are equal.This is further reflected in Fig. 6 where the regressed models perform noticeably better than those developed in this work.
Comparing both Eriksen et al. 's and Selam et al. 's models, in terms of the %AAD, it would seem that the former is moreaccurate.However, visually inspecting Fig. 6, 7 a and 7 c, the situation appears to be inverted.The reason for this again relates to the underlying model used for the RSP; the conditions at which the %AADs are determined approach the point at which Zuo and Fürst [63] 's correlation approaches zero and, as shown in Fig. 5 a, the impact of this is felt even at temperatures below this point.We can see this visually in figures 7 b and 7 d where Selam's model observes significant deviations from experimental data.
It is apparent from Table 10 that, again, the RSP model has the most-significant impact on the performance of a given electrolyte model in predicting these properties, where the MIAC appears to be more-sensitive to this than the osmotic coefficient.This is also reflected in Fig. 8 where we observe drastically different predictions between the different models.Interestingly, the electrolytedependent models (5 and 6b) appear to shift the MIAC and osmotic coefficient upwards compared to solvent-dependent models (1, 4 and 6a).Whilst all models observe an initial decrease in MIAC, the approach which used a constant permittivity (1) did not observe the subsequent increase, resulting in the minima that is typically expected in this system.As we have shown in Fig. 6 with Selam et al. 's model, a composition-independent permittivity model does not necessarily imply that these cannot predict increasing osmotic coefficients.However, this increase is likely due to the contributions from the SAFT-VR Mie equation and not the Born and electrostatic contributions.This contribution in turn depends on the binary interaction parameter between water and ions, further emphasising the importance of this parameter.
Table 10 also seems to suggest that the performance of a given model does not depend significantly on the electrostatic interaction model used.This is further reflected in Fig. 6 where, for the same RSP model, the MSA model predicts osmotic coefficients slightly smaller than the DH model, with this gap between the models increasing with the molality of the electrolyte.
Finally, all models struggled to capture the pressure and temperature dependencies of these properties.All predicted MIAC and osmotic coefficients are visually independent of pressure in Fig. 7 a and 7 c.We note that if the model does not predict an accurate MIAC or osmotic coefficient at atmospheric conditions, this remains a practically systematic error at all other pressures.In terms of the temperature dependence, we do observe the repercussions of the root present in Zuo and Fürst [63] 's correlation as Selam et al. 's model begins to diverge significantly from experimental values, accounting for the poor %AADs in Table 10 .Nevertheless, with the exception of models using Maribo-Mogensen et al. [75] 's permittivity model (model 6b), all other models did not predict the presence of a maxima in the MIAC and osmotic coefficient.Only models which use Maribo-Mogensen et al. [75] 's permittivity model including electrolyte association predict the presence of a maximum at higher temperatures.Given that neither Eriksen  In general, it is likely that the poor performance of these models at conditions other than room temperature and pressure is due to the fact that the model parameters are optimised without using data from these more-extreme conditions.This perhaps highlights the importance of using data from a wide range of conditions to train the model parameters.[107] of the saturation pressure of an aqueous solution of sodium chloride predicted using different electrolyte SAFT-VR Mie models and RSP models.Temperature range of 294 K to 367 K and molality range of 1.0 mol kg −1 to 5.0 mol kg −1 .

Solvent saturation pressure
A property that is generally of greater interest for electrolyte systems is the influence of the electrolyte on the saturation pressure of the solvent.The properties discussed in the previous section only examined a single phase; here, we examine the influence of electrolytes on a system at equilibrium with two phases, thus involving derivatives in both volume and composition.As mentioned in section 2.5 , we have restricted the electrolyte to only be present in the liquid-phase.This does not affect models which use density-dependent RSP models as the ions should be driven towards the high-density phase.%AAD are given in Table 11 , with the composition-dependence of the solvent saturation pressure predicted using some of models developed illustrated in Fig. 9 .
Indeed, as it is for most properties, the RSP model appears to have the most-significant impact on the performance of a given electrolyte model.Models that use either the DH or MSA term have similar predictions, with the difference between them increasing with molality.However, in comparison to other properties examined in this study, the extent of this influence is quite limited.All %AADs in Table 11 are smaller than what is observed for other properties.For similar reason to the poor performance of the second derivative properties in section 3.3 , here, the water model was regressed using saturation pressure data over a range of temperatures.Furthermore, the addition of any solute with unfavourable interactions with water (as is the case with ions) will lead to a reduction in saturation pressure.As a result, even without the influence of the electrolyte terms, we will observe the correct trend in saturation pressure against molality of the electrolyte.This is indeed what we observe in Fig. 9 when using a constant RSP with the Debye-Hückel model; in this case, the only terms having any influence on the saturation pressure are those relating to the SAFT-VR Mie equation.
Even when including the solvent composition dependence in our RSP model, the impact on the saturation pressure is limited.This is observed in both Table 11 , where we observe only a 1% decrease in %AAD, and Fig. 9 where, even at higher molalities, only a 0.004 MPa difference is observed between models using a constant and solvent-composition dependent RSP.Interestingly, the inclusion of the electrolyte-composition dependence through Maribo-Mogensen et al. [75] 's model (6b), a large underestimate of the solvent saturation pressure is observed, highlighting the significant impact of including this dependence in our RSP model.
Between the electrostatic interaction models, it seems that, in most cases, they have a limited influence on the saturation pressure.However, in models that depend on the electrolyte composition, we observe larger differences in %AADs.We can see the reason for this visually from Fig. 9 where differences in the electrostatic interaction models increase at higher molalities (similar to what is observed in Fig. 6 ).
From the previous discussion, it is apparent that the SAFT-VR Mie contributions have the largest impact on the saturation pressure.This is further highlighted by the improved performance of both Eriksen et al. 's and Selam et al. 's models which regressed the unlike potential depth parameter between water and the ions.Surprisingly, although Eriksen et al. 's model used saturation pressure data to train its model parameters, it is Selam et al. 's model which observes the lower %AADs, although the latter model allowed for two adjustable parameters whereas the former only regressed the unlike potential depth parameter.

Solubility of carbon dioxide
With the growing interest in carbon-capture processes, it is typically desired to develop extensible electrolyte models to predict the solubility of carbon dioxide [1,43] .However, a barrier to this is the fact that the only data available to regress the binary interaction parameter between carbon dioxide and ions is the solubility data which we hope to predict.As of the time of writing this study, no unlike interaction parameters between carbon dioxide and the sodium and chloride ions have been developed for either Eriksen et al. 's or Selam et al. 's model.Furthermore, for the water model used in Selam et al. 's work, no unlike interaction parameters between water and carbon dioxide have been developed.As a result, neither of these models can be used in predicting the solubility of carbon dioxide in electrolyte solutions.
However, it is possible to use the combining rule given in Eq. ( 63) to obtain the unlike potential depth parameter between carbon dioxide and ions a priori .As carbon dioxide only has a quadrupole, many of the terms in Eq. ( 63) are eliminated (giving CO 2 −Na = 80 .250K and CO 2 −Cl = 162 .85K when r = 1 ).Unlike σ i j and λ k,i j are obtained using the standard combining rules.
An additional limitation is the RSP models.Naturally, we can use a constant RSP or Zuo and Fürst [63] 's correlation, neglecting the influence of carbon dioxide.We can also use Schrecken- berg et al. [43] 's and Maribo-Mogensen et al. [69] 's solvent correlations as these can easily account for the influence of carbon dioxide.The latter is particularly interesting as, within the SAFT-VR Mie framework, carbon dioxide can associate with water, influencing the latter's association fraction, which is then used to predict the RSP.If we wish to extend to more-complex models, we find that Michelsen and Mollerup [68] and Maribo-Mogensen et al. [75] electrolyte correlations cannot readily be extended to carbon dioxide-containing systems.We will also not be investigating the use of Valiskó and Boda [64] 's correlation as, not only was it not regressed to systems containing carbon dioxide, but the conditions we will be investigating are beyond those in which the correlation was regressed.
Nevertheless, in Fig. 10 , we observe the predicted carbon dioxide solubility using the parameters developed in this study when implemented in different electrolyte models.The agreement with experimental data is quite limited, likely due to the un-optimised parameters.However, we do note that, even in the electrolytedevoid system, the agreement with experimental data is poor.
Between electrostatic interaction models, there is little difference, visually, between figures 10 a and 10 b.Once again, we can see that it is the RSP model that has the largest influence on the predicted property.Moreover, in this case, the influence is quite drastic as, between using a constant and solvent-dependent RSP, the phenomenological behaviour of the system changes from the electrolyte having a salting-in to a salting-out effect.This emphasises the importance of selecting a reliable RSP model prior to regressing model parameters.
It is difficult to assess whether including the electrolyte composition dependence would improve the model predictions.Nevertheless, increasing the composition of the electrolyte should subsequently reduce the RSP of the medium, reducing the stability of the ions within the solvent and, thus, increasing the bubble pressure.This should improve the agreement with experimental data.As such, a electrolyte-composition dependent RSP model that can be extended to carbon dioxide would be of interest for those wishing to obtain better agreement with experimental data.However, a similar effect could be achieved by adjusting the binary interaction parameters between carbon dioxide and the ions to make such interaction more unfavourable.As mentioned earlier, this would require the use of the same data we are trying to predict in order to regress the parameters.
Overall, it is apparent that the RSP model has a very significant influence on the solubility of carbon dioxide within an electrolyte solution and careful considerations should be made when developing a model to predict such properties.

Discussion
In this section, we take the opportunity to discuss in greater detail two aspects of this study.Firstly, we examine the contributions of each term in the electrolyte equation of state to the chemical potential of each species in the NaCl( aq ) system when using different ionic and RSP models.Subsequently, we examine the effectiveness of the proposed combining rule for electrolyte systems other than NaCl( aq ).

Contributions within an electrolyte equation of state
Throughout section 3 , we have examined the overall properties of an electrolyte system and, particularly with respect to the osmotic and mean ionic activity coefficient, we make reference to the relative contribution of each of the terms in the electrolyte equations of state examined in this work.In this section, we analyse these contributions in greater detail, and how these change when different RSP models are used.
As we have established in section 3 , it is the properties related to the chemical potential (or activity coefficient) of each species which are most-affected by these different models.As such, we will focus solely on the contributions of each term to the activity coefficient ( ln γ = (μ res ,i − μ * res ,i ) / (k B T ) ) for each species in the NaCl( aq ) system.However, for different combinations of ionic and RSP models, we obtain different densities, which, in turn, result in different values of the activity.Visualising this would require many curves on a single figure, each of which are likely going to be similar given the nominal change in density with salt molality is quite minor (as shown in Fig. 4 ).As a result, similar to the work by Maribo-Mogensen et al. [41] , we will instead examine the NaCl( aq ) system with increasing salt molality at 298 K and the density of pure water at atmospheric pressure determined using SAFT-VR Mie.The results are shown in Fig. 11 .
In general, the MSA and DH models, for the same RSP model, predict similar activities (the extent of this similarity does depend slightly on the RSP model); this is indeed what we have observed throughout section 3 .Interestingly, for the ion-related properties, the difference in electrostatic interaction terms when implementing a solvent-composition dependent RSP model compared to a constant RSP is quite small.In both cases, the contribution to the activity appears to saturate with increasing molality.This is likely to be related to the inverse screening length present in both electrostatic interaction models given in Eq. ( 14) and (19) .With increasing molality of salt, we observe reductions in the screening length as a result of the presence of more ions and, in the However, when the electrolyte-composition dependence is included within the RSP model, we observe a significant decrease in RSP, resulting in further decreases in the screening length from which an increase in the magnitude of its contribution to the activity of the ions is expected.This is likely to be magnified by the fact that, in both the DH and MSA models, there is a pre-factor of −1 r .This would imply greater stability, which is not what we would typically expect with increasing molality as, with a higher packing of ions, we would expect unlike interactions to begin dominating and, eventually, precipitation or phase separation.
On the other hand, the Born solvation term exhibits the opposite trend compared to the electrostatic interaction terms for the ionic species.With a constant and solvent-composition-dependent RSP model, the Born term has an almost insignificant contribution to the overall property.This is perhaps unsurprising as the contribution to the activity of the Born term is primarily related to the pre-factor 1 r − 1 .In the case of a constant and solventcomposition-dependent RSP, the change in this will be very small (or zero in the case of the constant RSP).However, with an electrolyte-composition-dependent RSP model, this pre-factor will increase, resulting in the behaviour observed in figures 11 a and 11 b.Now examining the solvent in Fig. 11 c, similar behaviour to the ions is observed.When a constant RSP is used, neither the Born solvation, DH or MSA term contribute to the activity; this is unsurprising as none of these terms carry an explicit solvent dependence beyond the RSP.When a solvent-composition-dependence is introduced to the RSP, in contrast to the ions, a significant change in the contribution of both the Born and electrostatic interaction terms is observed for the activity of the solvent.The sign of these contributions has also been reversed, compared to the ions, where the ionic and Born terms now have positive and negative contribution to ln γ , respectively.As it is with the ions, this is further amplified when using an electrolyte-composition-dependent RSP.This behaviour is to be expected as the solvation of ions allows for the formation of water-ion clusters, resulting in favourable interactions which reduce the activity.On the other hand, by accounting for the electrostatic interactions through the DH or MSA terms, we see a positive contribution to ln γ of the solvent.This arises from the reduction in RSP through which the solvent dependence arises in the electrostatic interaction terms; this reduction then arises from, at least within the RSP models considered [69,75] , the re-orientation of the solvent dipoles around the center ion (other reasons may include kinetic depolarisation [25] but these effects are not accounted for in the models considered).
For both ions and the solvent, the SAFT-VR Mie contribution appears quite small when examining figures 11 a to 11 c.For the cation and the solvent, this contribution is positive, which is to be expected considering that the hard-sphere, dispersive and associative interactions between these species is likely to be unfavourable.Surprisingly, this is not the case for the anion which has a slightly negative contribution.Whilst this could be as a result of poorly-estimated parameters, Eriksen et al. 's unlike interaction parameter between chloride and water is only approximately 10 K greater than the one obtained by the proposed combining rule.This may be due to the similar segment size of both water and chloride allowing for more-favourable interactions.Nevertheless, these favourable interactions are slight, relative to the other contributions to the activity.
Examining figures 11 a to 11 c may give us the false impression that the SAFT-VR Mie (or dispersive+associative) contributions are not of particular importance in such systems.However, as noted previously, for all species, the Born and electrostatic interaction terms are of opposite sign; when these are added to-Fig.12. Predicted mean ionic activity coefficient for various electrolytes using either the DH or MSA term, coupled with either a constant permitivity, the Schreckenberg or Maribo-Mogensen (Solvent, S and Electrolyte, E) models.Empty squares represent experimental data [113,[114][115][116][117][118] .
gether in figures 11 d to 11 f, their magnitude becomes comparable to that of the SAFT-VR Mie contribution.Whilst this does illustrate that the SAFT-VR Mie contribution should not be neglected, it also highlights the balance that is struck between the Born and electrostatic interaction terms.The importance of such balances have been discussed extensively by others in the field of electrolytes [24,30,[110][111][112] .However, we can now observe the role of the RSP in this balance.For the ions, in the case of a constant or solvent-composition-dependent RSP, it is clear that the electrostatic interaction terms are dominating and, without the SAFT-VR Mie contribution, we do not observe the expected increase in activity at higher molalities.However, once an electrolyte-compositiondependent RSP model is used, this balance is reversed where, although the electrostatic interaction terms initially dominate, at higher molalities, the Born solvation contribution dominates, resulting in the expected increase in activity.In the case of the solvent, the Born solvation term is clearly always dominating; to a certain extent, this is to be expected as the solvation effects are likely to be more significant than the charge screening effects for the solvent.It would be interesting to consider, in a future study, how the contribution from the NP-MSA model, which already accounts for the 'Born' energy [58] , compares to the MSA+Born contributions paired with different RSP models shown in Fig. 11 .This might give some indication as to which implementation of the RSP is most 'correct' and whether the representation of the ion-dipole interactions is captured adequately through the dispersion terms in the SAFT-VR Mie equation.
Overall, many of the observations made in section 3 can be explained by the behaviour observed in Fig. 11 .It is further emphasised that the difference between the DH and MSA models is quite negligible and it is the RSP model used that is of greater importance.It is also important to note that the role of the SAFT-VR Mie (and other dispersive and associative models) and the Born solvation term cannot be understated in electrolyte systems.In the case of the latter, the RSP model chosen can drastically change its overall contribution; this is quite surprising as other electrolyte models often neglect its contribution [26,42,74] although these do not use an electrolyte-composition-dependent RSP model.

Efficacy of proposed combining rule
Within this study, we have introduced a novel combining rule between ions and neutral species given in Eq. ( 63) .However, up until this point, we have only examined the efficacy of this equation for the NaCl( aq ) system where the results have been mixed (see section 3 ).As a result, we now extend and examine, in brief, the effectiveness of Eq. ( 63) in modelling other electrolyte systems.In order to cover a range of electrolytes (1-1, 1-2, 1-3, 2-1 and 2-3 types), we have limited ourselves to examining the performance of this equation when predicting the mean ionic activity coefficient.Still using r = 32 , we can obtain parameters for a variety of electrolytes which are given in Table 6 .The results are given in Fig. 12 where we still investigate the effects of using different electrostatic and RSP models.
We can begin to see issues with using the same RSP for different ions in Eq. ( 63) as, for the same stoichiometric coefficients, the predicted MIACs in figures 12 a to 12 h and 12 m to 12 p are almost indistinguishable.It is also clear that the role of the RSP model is still significant as, in some cases, it greatly improves the predictions of the electrolyte model (examples include MgCl 2 in Fig. 12 f with the DH model and AlCl 3 in Fig. 12 k with the MSA model).This does illustrate the need to tune r to the ion or electrolyte being modelled or regress the model parameters.
Another interesting observation is that this combining rule appears to struggle most with simpler salts (such as 1-1 salts).However, with more-exotic salts where experimental data are scarce (and where we would be more-likely to need a combining rule), the parameters obtained from Eq. ( 63) provide excellent agreement with experimental data (particularly for AlCl 3 , Al 2 (SO 4 ) 3 and La 2 (SO 4 ) 3 ).The best agreement is achieved when using solventcomposition-dependent RSP models such as Schreckenberg et al. [43] 's or Maribo-Mogensen et al. [69] 's.The exception to this is the Na 2 SO 4 electrolyte where, regardless of which electrolyte model is used, the model predictions are in strong disagreement with the experimental data.Whilst we could blame a poor guess of r , this disagreement could also arise from the omitted quadrupolar or associative interactions between water and the sulfate ion, which would be more significant than in the case of other ions.The reason this is not observed in the predictions for Al 2 (SO 4 ) 3 and La 2 (SO 4 ) 3 is that we only investigate low molalities for these electrolytes.
Overall, whilst the proposed combining rule shows promise, particularly with electrolyte systems where experimental data might be scarce, it is apparent that the RSP used has a large impact on the efficacy of the combining rule.It may be possible to leave r as an ion-specific or salt-specific adjustable parameter, where the latter may reduce the number of parameters needed to regress.Certainly the parameterisation of r warrants further studies beyond the scope of this work.

Concluding remarks
The problem of modelling the thermodynamic behaviour of electrolyte systems has sparked a great-many debates in the scientific community, with no unanimous agreement as to the 'correct' modelling strategy.This has resulted in multiple electrolyte equations of state used to model such behaviour where, even from a single underlying framework, many different variants can arise.Within this study, in the context of the SAFT-VR Mie framework, we examined and compared the impact of two aspects where such approaches typically differ: the model characterising the electrostatic interactions (Debye-Hückel and Mean-Spherical Approximation theory) and the underlying relative static permittivity model (constant, temperature-, density-, solvent compositionand electrolyte composition-dependent).The SAFT-VR Mie framework was selected primarily because of the existence of two approaches to model electrolyte systems which differ in both of the aforementioned aspects: the SAFT-VRE Mie [31] and eSAFT-VR Mie [32] equations.This allowed us to also consider the importance of parameter estimation in the modelling of electrolyte properties.
In order to allow for a fair comparison between different combinations of ionic and RSP models, the same underlying model parameters for a given electrolyte are used.Whilst most of the parameters are obtained either from experimental data or preexisting combining rules, to obtain the unlike interaction parameter between ions and neutral species, a novel combining rule was developed which account for the ion-dipole interactions.This combining rule has shown to provide good estimates of this parameter, particularly for electrolyte systems where little data are available.However, further analysis will be required to evaluate its efficacy in describing electrolyte systems globally; particular considerations will need to be made with respect to the relative static permittivity used to characterise the unlike interactions.
Using the NaCl( aq ) system as a case study, the impact of the ionic and RSP model on a variety of electrolyte properties was explored.The relative static permittivity, a property that is both predicted and used by the thermodynamic model, was found to be unaffected by the electrostatic interaction model used and the model parameters, even in cases where the RSP model is dependent on variables determined by the model parameters.The Maribo-Mogensen et al. [75] model for electrolyte systems, which carries a temperature, volume and composition dependence, proved to give the best overall agreement with experimental data.
The liquid densities proved to be effectively independent of the electrostatic interaction model used, with the relative static permittivity model having a slightly larger impact.However, it is the model parameters that had the biggest impact.Given that eSAFT-VR Mie resulted in the best agreement with experimental data, it could be hinted that the size of the ion has the larger impact on the density relative to the other parameters, considering the latter was treated as an adjustable parameter within the eSAFT-VR Mie framework.Similarly, second derivative properties are more so dependent on the relative static permittivity model used than the electrostatic interaction model, where solvent-compositiondependent RSP models yielded better performance.They are also relatively unaffected by the model parameters as even the regressed parameters observed a similar performance as the nonregressed ones.We also note that the performance of these properties is, for the most part, dependent on the accuracy of the solvent model used to predict such properties (inaccuracies in the solvent model lead to practically-systematic errors in the electrolyte system).
The osmotic and mean ionic activity coefficient are two properties which are affected by both the ionic and RSP model (more so the latter, however).It is unclear as to which RSP model is more 'correct'.We note that more-accurate estimates are obtained when using a solvent-composition-dependent model over an electrolytecomposition model, although, based on the results obtained by the Eriksen et al. [31] and Selam et al. [32] equations, it is clear that the model parameters used are the dominating factor in the resulting model performance.Nevertheless, in-spite of this, many of the models struggled to capture the temperature dependence of these properties, where only those models which used an electrolyte-composition-dependent RSP observed a maximum in the MIAC.This may indicate that such a RSP model would be more-appropriate, although, it may also be worthwhile examining the performance when above-ambient temperature activity coefficient data are used to regress the model parameters.
The solvent saturation pressure was found to be primarily dependent on the underlying SAFT-VR Mie model and, as a result, the model parameters, although the RSP model used also had a noticeable impact on this property.Like the second-derivative properties, the solvent saturation pressure is largerly dependent on the underlying solvent model.The final property examined was the carbon dioxide solubility in an electrolyte solution and, although there were limitations on which RSP models could be examined, it is clear that the RSP model can have an impact on the predicted thermodynamic behaviour as significant as changing the effect of the salt from salting-in to salting out.
Finally, we examined the relative contribution of each of terms present in an electrolyte equation of state to the activity of each species in an electrolyte system.Whilst much of the analysis corroborated what had been observed in the results (such as the small difference between electrostatic interaction terms), it also highlighted the importance of the Born-solvation term within such models and how it is drastically affected by the RSP model used.The balance between the different terms is also dependent on the chosen RSP model as, in the case of the ionic species, this can shift the balance between a dominating electrostatic interaction or Born solvation term.
Overall, although it is not the purpose of this study to identify the 'best' combination of electrostatic interaction and RSP models to develop an electrolyte model, we can draw some conclusions as to which approach might be more suitable.Firstly, it is clear that, between the Debye-Hückel and Mean-Spherical Approximation terms, the impact on the overall model is insignificant enough to make these almost interchangeable (especially using regressed model parameters).Secondly, although it is unclear as to whether the relative static permittivity should carry a composition dependence, at the very least, it should not be treated as a constant (special care should also be taken to ensure that the model used can extrapolate to more-extreme conditions without incurring unphysical values).Thirdly, if the RSP carries a composition-dependence, the Born-solvation term should not be neglected.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. i j between water and various ions predicted using the proposed combining rule for various ions compared to those estimated by Eriksen et al. and Selam et al. .
i j is now comparable to those obtained by Eriksen et al. and Selam et al. .More importantly, the trend between the molar mass of ions and i j is now very similar between those obtained from Eq. (63) and Selam et al. .Whilst the trend for cations is similar between Eq. (63) and Eriksen et al. , for anions, the trend is the opposite of what is obtained from Eq. (63) (particularly when looking at the fluoride ion).This may be due to the strategy used by Eriksen et al. (simultaneous optimisation) compared to the more systematic approach used by Selam et al. .

r)
dependence to the RSP model, the estimates of the models used in Fig. 4 increase to the point where the estimates of the DH model coupled with Maribo-Mogensen et al. [75] 's model (model 6b) in better agreement with experimental data than those made by Eriksen et al. 's model with regressed parameters.When compared to the estimates from Selam et al. 's model, which allowed σ ii to be an adjustable parameter, these estimates are far more-accurate than those made by any other model considered in Table

Fig. 5 .
Fig. 5. Predicted isobaric heat capacity ( figures 5 a and 5 b) and speed of sound ( figures 5 c and 5 d) using the SAFT-VR Mie+Born+DH equation of state.Empty squares represent experimental data [98,99] .

Fig. 6 .
Fig. 6.Predicted osmotic ( Fig. 6 a) and mean ionic ( Fig. 6 b) coefficients for NaCl( aq ) at p= 0.101 MPa and T = 298.15K using the SAFT-VR Mie + Born + DH (with adjusted or unadjusted σ ii ) or MSA coupled with the Maribo-Mogensen (solvent and electrolyte) correlation.Estimates from the Eriksen et al. and Selam et al. models have also been included.Empty squares represent experimental data [105] .

Fig. 7 .
Fig. 7. Predicted osmotic coefficient and MIAC for NaCl( aq ) at m = 6 mol kg −1 and T = 298.15K ( figures 7 a and 7 c) or p= 100 MPa ( figures 7 b and 7 d) using the SAFT-VR Mie + Born + DH or MSA coupled with the Maribo-Mogensen (electrolyte), Schreckenberg correlation or constant permittivity.Estimates from the Eriksen et al. and Selam et al. models have also been included.Empty squares represent experimental data [105] .

11 .
Chemical potential of water ( figures 11 c and 11 f) and the sodium ( figures 11 a and 11 d) and chloride ( figures 11 b and 11 e) ions in an aqueous sodium chloride mixture at a temperature of 298 K and molar volume of 0.018 dm 3 mol −1 (molar volume of pure water determined by SAFT-VR Mie at atmospheric pressure and temperature of 298 K). case of solvent-composition-dependent RSP model, decrease in the RSP of the medium.As more ions are introduced to the system, the screening length begins decreasing asymptotically, resulting in what is observed for a constant RSP in figures 11 a and 11 b.Even in the case of a solvent-composition-dependent RSP model, although a decrease in RSP also decreases the screening length [41] , as shown in Fig. 3 c, this decrease is quite negligible and accounts for the small differences observed in figures 11 a and 11 b.

Table 1
[47]ational temperatures and degeneracies of water and carbon dioxide obtained from Walker and Haslam[47].

Table 2 SAFT
[54,56] Parameters for the non-ionic species of interest in this study obtained from literature[54,56].

Table 4
Parameters for Eq.(29) corresponding to the various solvents of interest.

Table 6
i j / K parameter between water and various ions obtained from Eriksen et al. ,Selam et al. , Eq. ( [63]sen et al.to define phase equilibria.However, Selam et al. choose instead to constrain charged species such that they are not present in the gas phase, ignoring Eq. (76) .It is likely that the different approaches were used because of each model's chosen RSP model.By using Schreckenberg et al.[43]'s correlation within Eriksen et al. 's model, the difference in RSP between the two phases drives the ions into the liquid phase.In Selam et al. 's model, as Zuo and Fürst[63]'s correlation is used, which does not depend on the density of the phase, this driving force is not present which would result in a greater composition of ions in the gas phase.

Table 7
%AAD from experimental data

Table 8
%AAD from experimental data

Table 10
[105]from experimental data[105]of the Osmotic and Mean Ionic Activity coefficient of an aqueous solution of sodium chloride predicted using different electrolyte SAFT-VR Mie models and RSP models.Temperature range of 273.15K to 573.15 K, pressure of 0.101 MPa to 100 MPa and molality range of 0.0 mol kg −1 to 6.0 mol kg −1 .

Table 11
%AAD from experimental data