Publication Cover
Molecular Physics
An International Journal at the Interface Between Chemistry and Physics
Volume 114, 2016 - Issue 18: Thermodynamics 2015 Conference
2,761
Views
43
CrossRef citations to date
0
Altmetric
Articles

Development of intermolecular potential models for electrolyte solutions using an electrolyte SAFT-VR Mie equation of state

, , , , &
Pages 2724-2749 | Received 02 Mar 2016, Accepted 08 Sep 2016, Published online: 25 Oct 2016

ABSTRACT

We present a theoretical framework and parameterisation of intermolecular potentials for aqueous electrolyte solutions using the statistical associating fluid theory based on the Mie interaction potential (SAFT-VR Mie), coupled with the primitive, non-restricted mean-spherical approximation (MSA) for electrolytes. In common with other SAFT approaches, water is modelled as a spherical molecule with four off-centre association sites to represent the hydrogen-bonding interactions; the repulsive and dispersive interactions between the molecular cores are represented with a potential of the Mie (generalised Lennard-Jones) form. The ionic species are modelled as fully dissociated, and each ion is treated as spherical: Coulombic ion–ion interactions are included at the centre of a Mie core; the ion–water interactions are also modelled with a Mie potential without an explicit treatment of ion–dipole interaction. A Born contribution to the Helmholtz free energy of the system is included to account for the process of charging the ions in the aqueous dielectric medium. The parameterisation of the ion potential models is simplified by representing the ion–ion dispersive interaction energies with a modified version of the London theory for the unlike attractions. By combining the Shannon estimates of the size of the ionic species with the Born cavity size reported by Rashin and Honig, the parameterisation of the model is reduced to the determination of a single ion–solvent attractive interaction parameter. The resulting SAFT-VRE Mie parameter sets allow one to accurately reproduce the densities, vapour pressures, and osmotic coefficients for a broad variety of aqueous electrolyte solutions; the activity coefficients of the ions, which are not used in the parameterisation of the models, are also found to be in good agreement with the experimental data. The models are shown to be reliable beyond the molality range considered during parameter estimation. The inclusion of the Born free-energy contribution, together with appropriate estimates for the size of the ionic cavity, allows for accurate predictions of the Gibbs free energy of solvation of the ionic species considered. The solubility limits are also predicted for a number of salts; in cases where reliable reference data are available the predictions are in good agreement with experiment.

1. Introduction

Aqueous electrolyte solutions are liquid mixtures containing charged species, typically stemming from the dissociation equilibrium of the constituent ions of salts in water. The ubiquitous nature of electrolytes makes them relevant in many scientific and industrial applications. Examples of applications requiring a thermodynamic description of electrolyte solutions include the sequestration of carbon dioxide, the clean-up of process chemicals from nuclear waste, hydrate formation, scaling in pipelines, advanced process chemistry, corrosion processes, and flue-gas cleaning. Electrolyte solutions are also prevalent in geochemical and biological environments and the development of molecular models to describe their properties remains of great interest.

The presence of charged species (ionic components) in solution drastically alters the bulk thermodynamic properties of the mixture, usually leading to a lowering of the vapour pressure (with a corresponding increase in the boiling temperature), an increase in density, and changes in the viscosity and surface tension of the system. The theoretical description of thermodynamic properties of electrolyte solutions has proved a challenging task, and the body of work in this area has been less extensive than for non-electrolyte systems. The development of a successful theory of electrolyte solutions relies on both the availability of accurate inter-particle potentials for all of the species in the system and the ability to resolve the statistical–mechanical relations that lead to the bulk properties of interest. Much of the difficulty can be ascribed to the range of the electrostatic forces, and the corresponding paucity of analytic theories that can be used to describe the relevant many-body effects. Available theories predominantly incorporate the Coulombic nature of the ionic interactions with a much simplified representation of the neutral solvent molecules, which are often polar and, as in the case of water, can exhibit hydrogen-bonding interactions. The intricacy of the solvent–solvent, solvent–solute and solute–solute interactions and their corresponding interplay makes this a challenging area of research.

A so-called primitive formulation is proposed in the classical approach of Debye and Hückel (DH) [Citation1], whereby the solvent is represented as a continuum dielectric medium and the Poisson–Boltzmann equation [Citation2,Citation3] is solved in linear form. Waisman and Lebowitz [Citation4] and Blum and Hoye [Citation5,Citation6] later solved an equivalent primitive model using integral equation theory within the mean-spherical approximation (MSA). A formulation of the MSA allowing for an explicit representation of the solvent was then developed by Wei and Blum [Citation7], among others. In the late 1990s, molecular-based equations of state (EOSs), specifically the statistical associating fluid theory (SAFT) [Citation8,Citation9] and the cubic plus association (CPA) EOS [Citation10], were combined with the classical approaches for charged species [Citation1,Citation5,Citation6] and shown to provide an accurate description of the properties of strong electrolyte solutions [Citation11–13], extending the range of applicability to more-concentrated solutions. The increasing accuracy in the description of the properties of the solvent and fidelity of the representation of the solvent–ion interactions has been exploited in subsequent work [Citation14–22].

For applications where the partitioning of the charged species in different phases is of interest, further consideration of the electrostatic contribution that arises from the interactions of the ionic species with the medium is required. In the case of non-primitive models, where the polarity of the solvent is treated explicitly by including dipolar interactions, such a contribution arises naturally [Citation23,Citation24]. The solvent is treated explicitly only in terms of its non-electrostatic nature for most studies with the SAFT or cubic EOS electrolyte framework and, as a consequence, an additional contribution is necessary for an accurate representation of the thermodynamics (in particular, the chemical potential) of the ionic species due to the dielectric nature of the real medium. A simple treatment to account for the introduction of a charged species in a dielectric medium was provided by Born [Citation25]. The Born contribution has now been incorporated within cubic and CPA EOSs [Citation26,Citation27], and, more recently, within the ePPC-SAFT [Citation20] and SAFT-VRE [Citation21] approaches. A non-primitive treatment that incorporates a polar solvent within a SAFT formalism has also been presented [Citation17,Citation28–30]. These studies highlight the importance of accurately predicting the dielectric constant of the system, as well as the impact of taking into account ion-pairing phenomena. The work of Maribo-Mogensen and co-workers [Citation22,Citation31,Citation32] is also particularly relevant in this context, where the focus is on improving the description of the static permitivity, a key property of electrolyte solutions and which is commonly incorporated empirically from a knowledge of the experimental value.

The SAFT-VRE [Citation13] approach has recently been revised [Citation21] to acknowledge the specific importance of incorporating the Born [Citation25] contribution to the free energy of the system, thus accounting for the charging process of the ionic species in the solvent medium. A robust set of models for strong electrolytes was presented, providing an accurate description of a broad range of properties, including the Gibbs free energy of solvation, with the aid of the Born treatment. The model and theory proposed in [Citation21] are based on an underlying square–well interaction potential following the original SAFT-VR framework [Citation33,Citation34]. The SAFT-VR EOS has now been reformulated for the more versatile and realistic Mie [Citation35] potential form together with an improvement in the underlying statistical–mechanical theory, now taken to third order in the high-temperature perturbation expansion [Citation36]. The novel SAFT-VR Mie EOS has brought a new level of reliability in the development of an accurate platform for the thermodynamic properties of complex fluids and fluid mixtures [Citation36–38]. The approach has also been cast as a group-contribution methodology (SAFT-γ) [Citation39,Citation40], where the molecules are represented in terms of the underlying chemical functional groups, providing an enhanced predictive capability. Importantly, the improved theoretical framework provides models that reflect more accurately the physics of the systems they describe, reducing the reliance on re-adjusting the values of the intermolecular parameters to capture given thermodynamic properties of interest (although, of course, the estimation of the key parameters is still necessary). These advances have not yet been extended to include a description of electrolyte solutions. The purpose of our current work is to redress this: as will be shown, the introduction of a more rigourous physical description in the resulting SAFT-VRE Mie EOS allows for a reduction in the number of adjustable parameters, as compared with the previous SAFT-VRE formulation [Citation21]. In addition, we develop electrolyte model parameters to represent the thermodynamic properties of a broad range of aqueous solutions of strong electrolytes including monovalent and divalent ions, and assess the performance of the models by making appropriate comparisons with available experimental data.

Our paper is set out as follows: in the following section, we describe the underlying SAFT-VRE Mie theory; in Section 3, the technique for the development of our model intermolecular parameters is discussed in detail; in Section 4, we present the models developed, which comprise group-I metal cations, from lithium to rubidium, selected group-II (divalent) metallic cations, the hydronium cation, the halide anions from fluoride to iodide, and the hydroxide anion; our models and theory are applied in Section 5 for the description of thermodynamic properties, including the mean-ionic activity coefficients, the osmotic coefficients, the vapour pressure, the density, and the solvation energies of the ions for aqueous solutions containing these species; we present our conclusions in Section 6.

2. Methodology

2.1. Thermodynamic system

At a specified temperature and total volume, we consider fluid mixtures consisting of water (the solvent) and a salt which is fully dissociated into its constituent ions following an equilibrium reaction. For a salt MX, the equilibrium can be described as (1) where the salt dissociates into the metal ion M of valency Z+ and the counterion X of valency Z in the stoichiometry ν+: ν. The dissociation reaction is fully shifted towards the products, i.e. no neutral salt is present in the mixture. Here we consider only the dissociated ions present in the solvent as solutes.

2.2. SAFT-VRE Mie Helmholtz free energy

In broad strokes, the theory as outlined in the present paper follows the approach developed in previous work [Citation21], as applied to the SAFT-VR Mie formalism [Citation36] in conjunction with the Wertheim [Citation41–46] TPT1 treatment of the association in water [Citation37] (based on the Lennard-Jones form of the association kernel). We refer to the resulting theory as SAFT-VRE Mie. A short overview of the main expressions is given here for completeness; for further details, the reader is referred to the original publications.

The Helmholtz free energy A of the system is written as a sum of terms, arising from a perturbation approach, each of which corresponds to a specific physical contribution in addition to the usual ideal contribution represented as Aideal: (2) where the non-electrostatic and electrostatic contributions are treated separately. The residual non-electrostatic terms Amonomer, Achain, and Aassociation describe the common SAFT contributions: the effect of introducing monomeric spherical segments interacting through Mie potentials; the effect of bonding the monomeric spherical segments into molecular chains; and the effect of molecular association through short-ranged directional forces, such as those arising in hydrogen bonding. Collectively, these terms correspond to the non-electrolyte SAFT-VR Mie EOS [Citation36]. The electrostatic terms comprise Aion and ABorn, both represented within the primitive model formalism, where the solvent is described by a dielectric medium (it is important to note, however, that the solvent is treated explicitly in the non-electrostatic part of the free energy). The change in free energy associated with the process of charging the ions in the solvent is incorporated in ABorn following the Born [Citation25] model, and the effect of electrostatic interactions between charged species is described with Aion using the primitive non-restricted MSA [Citation5,Citation6] model.

2.2.1. Non-electrostatic contributions

The non-electrostatic contributions are modelled following the SAFT-VR Mie EOS [Citation36], in place of the SAFT-VR SW EOS applied in previous work [Citation21]. In the SAFT-VR Mie formalism, molecules are modelled as chains comprising mseg bonded spherical segments of diameter σ, interacting through a spherically symmetric Mie potential [Citation35,Citation36] characterised by a repulsive exponent λr, an attractive exponent λa, and a dispersive interaction energy ϵ. The Mie pair potential for the interaction between two unlike spherical segments labelled i and j can be expressed as (3) where rij is the distance between the centres of segments i and j. Additional off-centre, short-ranged square–well potential sites can be assigned to segments of particular species to mediate hydrogen-bonding associative interactions.

The ideal Helmholtz free-energy term of the mixture is given by [Citation47] (4) where kB is the Boltzmann constant, N is the number of molecules, nc is the number of species (recalling that in our work, we assume the salts to be fully dissociated so that the anion and cation are treated as separate species), ρi = Ni/V is the number density of component i at volume V, where Ni is the number of molecules of component i, xi = Ni/N is the mole fraction of component i, and Λ3i is the de Broglie volume for component i which includes all the relevant translational, vibrational and rotational contributions to the kinetic energy.

The monomer contribution of the Helmholtz free energy for the mixture of spherical segments interacting via Mie potentials is obtained as a third-order perturbation expansion following the high-temperature perturbation scheme of Barker and Henderson [Citation36,Citation48]: (5) where the reference hard-sphere term AHS is represented with the expression of Boublík and Mansoori [Citation49,Citation50] for mixtures of hard spheres. An analytical expression for the first-order mean-attractive energy A1 is obtained through the use of the mean value theorem and a mapping of packing fractions, following the original SAFT-VR [Citation33] approach, as applied for Mie fluids within a range of exponents of 5 ≤ λ ≤ 100 [Citation36]. The second-order fluctuation term A2 is calculated using the improved macroscopic compressibility approximation (MCA) of Zhang [Citation51] with the modification proposed by Paricaud [Citation52]. An empirical expression is used for the third-order perturbation term A3, with coefficients obtained to reproduce the critical point and fluid-phase equilibrium of Mie fluids [Citation36].

Although not relevant in the case of the molecules considered in our current work, for completeness and to provide a general framework, it is useful to mention briefly the contribution to the free energy due to the formation of chains of bonded Mie segments. This contribution, as in other SAFT approaches, is obtained through the TPT1 formalism of Wertheim [Citation46,Citation53], which is expressed in terms of the radial distribution function (RDF) of the monomer fluid, evaluated at contact distance of the segments. In the case of a mixture, the expression is given as (6) where is the RDF of the monomeric Mie fluid system at contact σii. Note that in our current work, the number of segments per chain is mseg, i = 1 for all of the species i considered and, as a result, Achain is equal to zero.

The association contribution of the Helmholtz free energy is also derived from the formalism presented by Wertheim [Citation41,Citation42], which can be written as [Citation45,Citation46] (7) for the case of a mixture of nc species, and where the second summation is over the types of sites of species i; nsites, i is the number of site types on a molecule and na, i is the number of sites of type a on component i. The variable Xa, i represents the fraction of molecules not bonded at a given site type a on species i and is obtained as the solution of a set of mass action equations for each bonding interaction: (8) In this expression, indices i and j represent the components, indices a and b represent the sites, and Δab,ij is the integrated association strength between site type a on molecule i and site type b on molecule j. The integrated association strength can be expressed [Citation37] as a product of the association kernel I, the Mayer function Fab of the bonding interaction of strength ϵHBab, ij, and the bonding volume Kab: (9) The association kernel I has been parameterised as a polynomial in the dimensionless temperature T* = kBT/ϵ, dimensionless density ρ* = Nσ3/V, and the repulsive exponent λr, as reported in [Citation37,Citation38]; note that in our current work, we employ the Lennard-Jones (12-6) form of the association kernel.

2.2.2. Electrostatic contributions

The Coulombic interactions pertaining to the charged species are represented following an unrestricted primitive model where each of the ionic species is considered to be spherical with a central point charge; the diameters of the ions are not restricted to be all identical. In this approach, the explicit representation of the polarity of the solvent molecules is replaced by a dielectric medium, of relative static permittivity D. Such a treatment of the solvent simplifies the inherently complex model for the Coulomb potential of the charged ions [Citation54], while retaining the main feature of electrolytes. The other ion–ion and ion–solvent non-electrostatic interactions (repulsion, dispersion and, if needed, hydrogen bonding) of the ionic species are taken into account using the SAFT-VR Mie expressions as outlined in Section 2.2.1.

The change in the free energy due to the insertion of all of the ions in the dielectric medium is accounted for by the free energy of solvation. In our work, the process of solvation is described with the Born approach [Citation25] which, though simplistic, captures the basic thermodynamic features. In the Born model, a spherical cavity is created, of a size suitable for each of the individual ions. A point charge of magnitude and sign commensurate with that particular ion is then embedded, and an expression for the free-energy change associated with this process is obtained through the Born cycle, which is illustrated in Figure : the system is first discharged to represent a corresponding system of non-interacting hard spheres in a vacuum; the dielectric medium is then switched on; finally, a non-interactive recharging of the ions is carried out. The premise of the Born model is that a spherical cavity of diameter is created in the dielectric medium for each ion i, independent of any others, leading to the following contribution to the Helmholtz free energy: (10) Here, the index i runs over all nion charged species, e = 1.602 × 10−19 C is the elementary charge, ε0 = 8.854 × 10−12 C2 J−1m−1 is the static permittivity of the medium in vacuum, D is the relative static permittivity of the medium, and Zi is the valency of ion i.

Figure 1. Schematic representation of the Born cycle. (1) The system consists of hard spheres and non-interacting charged hard spheres in vacuum; (2) the charged hard spheres are discharged in vacuum; (3) the system consists of (uncharged) hard spheres in a dielectric medium represented by the dotted background; (4) the system consists of hard spheres and non-interacting charged hard spheres in a dielectric medium.

Figure 1. Schematic representation of the Born cycle. (1) The system consists of hard spheres and non-interacting charged hard spheres in vacuum; (2) the charged hard spheres are discharged in vacuum; (3) the system consists of (uncharged) hard spheres in a dielectric medium represented by the dotted background; (4) the system consists of hard spheres and non-interacting charged hard spheres in a dielectric medium.

A particular challenge of implementing the Born model is the interpretation and parameterisation of the cavity diameter. The primitive model, in which the Born model is expressed, does not represent a convenient 1:1 mapping to the intermolecular potential parameters utilised in the explicit solvent representation of the system. It is common, however, to take the cavity diameter of the Born model as that represented by the diameter σii of the intermolecular potential of the ion itself.

Ion–ion interactions are represented with the Coulomb potential in a dielectric medium, i.e. (11) where rij is the centre–centre distance between charges i and j. In the case of the primitive model, one of two classical theories is typically used to resolve the corresponding residual free energy: either the Debye–Hückel theory [Citation1] or the mean spherical approximation (MSA) [Citation6,Citation55]. As shown by Maribo-Mogensen et al. [Citation56], both approaches lead to a similar representation of the macroscopic thermodynamic properties of electrolyte solutions; the material difference between the two is the parametrisation procedure for the corresponding models. In our current paper, the MSA theory is applied following the formulation of Blum [Citation5,Citation6], as laid out in previous work [Citation21]. The change in the Helmholtz free energy due to the electrostatic interactions between charged species within the MSA formalism can be expressed as (12) Here, UMSA is the MSA contribution to the internal energy U, and Γ is the screening length of the electrostatic forces. The MSA internal energy is given by (13) Δ represents the packing fraction of the ions as a function of their diameter σii: (14) The functions Pn and Ω are coupling parameters, where Pn couples to the charge of the ions, whereas Ω relates to the packing fraction of the ions. Both are functions of the ionic parameters, as well as the screening length of the ions: (15) (16) Finally, the screening length Γ is a function of the relative static permittivity and the effective charge Qi(Γ) of the ions, leading to an implicit formulation through Qi: (17) where the effective charge is related to the electric charge of the individual species and the Pn coupling parameter: (18)

The implicit nature of the MSA formulation requires an iterative procedure to establish Γ. In the present implementation, a successive substitution procedure is followed with an initial guess, Γ0, for the screening length obtained from the Debye–Hückel estimate for the screening length, (19) where κ is the inverse Debye–Hückel length. The addition of the Born term and the MSA term to the SAFT-VR Mie theory, detailed in Section 2.2.1, constitutes the SAFT-VRE Mie theory developed in our current work.

2.2.3. Auxiliary model: the relative static permittivity

The relative static permittivity D is calculated with a Harvey–Prausnitz model [Citation57] as outlined in a previous paper [Citation21]. This correlative model depends on the solvent composition, density, and temperature, and thereby exhibits an implicit dependence on the ion concentration. The value of D in the case of a pure solvent j is characterised by a volume parameter dV,j and by a temperature parameter dT,J. A combined static permittivity parameter dj is calculated for a given temperature following the relation (20) The relative static permittivity is then obtained as (21) where ρj = Nj/V is the number density of the solvent.

2.3. Thermodynamic properties of electrolyte solutions

The concentration of a given salt MX is often quantified in terms of the molality, denoted here by mMX, and defined as the molar amount (mol) of salt per kilogram of solvent j. The mole fraction of a given ion i can be calculated from the molality as (22) where νi is the stoichiometric coefficient of ion , and MWj denotes the molecular weight of the solvent in units of kg mol−1.

In addition to the common thermodynamic properties of equilibrium fluids, the activity coefficients γi of the ions and the osmotic coefficient Φ of the solution are key properties often reported for electrolyte solutions. The generic activity coefficient of component i, denoted here by γi, x, is related to the chemical potential μi of component i through the following expression: (23) the product ciγi, x yields the activity ai of component i. In Equation (Equation23), is a reference term, N is the composition vector, ci is a measure of concentration, and the subscript x denotes that γi, x is expressed in units commensurate with those of ci. Depending on the nature of the system in question, it may however be convenient to express the activity coefficient in one of a variety of ways. For the standard symmetrical activity coefficient, the concentration is measured in terms of the mole fraction xi of component i: (24) where the reference term is identical to that of the chemical potential of pure component i at the system temperature and pressure, i.e. , and the activity can be calculated as ai(T, p, N) = xiγi, x(T, p, N). As ions cannot be related meaningfully to the pure system, the appropriate reference for the activity coefficients is the rational asymmetric scale, for which the activity coefficient of ion i, γi, *, has a value of one in the limit of infinite dilution of ion i. In practice, the calculation of γi, * is carried out through the fugacity coefficient ϕi: (25) where N* is the composition vector at infinite dilution of ion i.

The fugacity coefficient is obtained from the residual chemical potential μRes = μ − μideal, taking into account the ideal state conversion for properties expressed at a given pressure, obtained from an EOS expressed with volume as the independent variable. In this case, the fugacity coefficient is obtained as (26) Here, Vp denotes the volume corresponding to the specified pressure and Z = pVp/(NkBT) denotes the thermodynamic compressibility factor.

It is the usual convention to report properties of ions in terms of the molality, whereby the activity coefficients are presented in the molal-based scale, γi, m. This scale employs a convention of a hypothetical ideal solution at unit molality such that ∑imi → 0 ⇒ γi, m → 1. The conversion between the rational asymmetric scale and the molal-based scale follows as [Citation58] (27) with a corresponding change in the reference term utilised in the calculation of the chemical potential based on γi, m, and where xj denotes the mole fraction of the solvent.

The electrolytic properties of the solution can be collated in the mean ionic activity coefficient (MIAC), γ±, m, calculated as an average of cationic and anionic contributions: (28) The contribution of the solvent j to the thermodynamics of the system can be characterised through the osmotic coefficient Φ. The osmotic coefficient provides a convenient re-scaling of the activity of the solvent, in order to accentuate its variation at low solute concentrations: (29) The properties of solutes and solvents are not independent, and must fulfil the Gibbs–Duhem relation [Citation59]. The MIAC and the osmotic coefficient Φ are related through the following relationship: (30)

Finally, Φ and γ±, m can be related to the change in the Gibbs energy of solvation ΔGsolv, i of ion i through the fugacity coefficient: (31) where ϕi is the fugacity coefficient of component i at infinite dilution in solvent j, pref is the pressure of the reference state for the change in the Gibbs free energy, and m° is the standard state molality of 1 mol kg− 1.

2.4. Phase equilibrium

The conditions which enforce phase equilibrium of electrolyte solutions are similar to those of non-electrolyte systems, with additional constraints related to the charge balances and the necessary charge neutrality of a given phase. First, the thermal and mechanical equilibrium conditions must be satisfied, i.e. (32) (33) Here, the index Nphases accounts for all of the phases, denoted as superscript Greek letters. In addition, a relationship in each equilibrium phase is required for the chemical potentials of each species in the mixture. As a consequence of treating electrolyte solutes as fully ionised in solution, the additional constraints required to characterise phase equilibrium will depend on the nature of the phases considered.

2.4.1. Fluid-phase equilibria

In the consideration of equilibrium between two (or more) fluid phases, equality of chemical potentials is required for each neutral species. Here it is assumed the solvent is pure and the only neutral species present so that: (34) while for each pair of charged species i and i′, a constant relative difference of chemical potentials across electro-neutral phases is satisfied [Citation60]: (35) (36)

The solution to the equilibrium conditions is obtained with a Levenberg–Marquardt [Citation61,Citation62] algorithm, allowing for the presence of salts in all fluid phases, including the gas phase. The volume dependence of the relative static permittivity and the inclusion of the Born free energy term in the underlying theory deliver a model where the ions naturally partition predominantly into the denser liquid phase, with only trace amounts in the gaseous phase, in agreement with experimental observation.

2.4.2. Solid–liquid phase equilibria

The description of phase equilibrium between solid and liquid phases containing electrolytes is linked to the chemical equilibrium governing the dissociation of the solvated electrolyte leading to the formation of charged species in solution. The solid phase consists of the pure unsolvated crystalline salt MX in equilibrium with a liquid phase saturated in salt. Assuming complete dissociation of the dissolved salt, the phase and chemical equilibria require that the chemical potential of the crystalline salt in the solid phase is equal to the sum of the chemical potentials of the solvated ions , with i =M,X: (37) where M represents the cation and X the anion. The molal composition vector of the saturated aqueous phase is represented by msat, and the molality of solvated ion i is related to the salt molality through (38) where mMXsat is the solubility limit of salt MX.

The chemical potential of the solid salt is obtained from (39) where μ0MX is the reference chemical potential of the pure solid salt and aMX is the activity of the pure salt, taken to be unity; the reference chemical potential of a pure compound is also its molar Gibbs free energy of formation .

The chemical potential of solvated ion i is expressed on a molality basis as (40) where the reference chemical potential μ°i(T, p, m°) of ion i refers to a hypothetical ideal solution of unit molality (m = 1 mol kg−1), and γi, m is the molal-based activity coefficient given by Equation (Equation27). The chemical potential of ion i at the reference state of unit molality corresponds to the molar Gibbs free energy of formation of 1 mol kg−1 of the solvated ion.

Using Equation (Equation38) and the expressions for the chemical potential of the solid salt and solvated ions, the solid–liquid equilibrium condition for fully dissociated salts given by Equation (Equation37) can be rewritten to obtain a solubility equation for the salt as (41) In this expression, Ksp is the solubility product of the salt given by (42)

3. Model development

We consider aqueous solutions of strong electrolytes focusing particularly on the halide salts of alkali metals and alkaline earth metals, as well as aqueous strong acids and bases. These solutions are modelled as ternary mixtures composed of water, anions, and cations, under the assumption of a fully dissociated solute. The electrolyte solutes are therefore modelled via the constituent monovalent and divalent atomic ions in the case of salts, while molecular ions are also considered for the acid and base solutions.

In our current work, the ions are modelled as spherical, consisting of a single segment mseg, i = 1) of diameter σii, and carrying a single point charge qi = Zie. All ions experience dispersion interactions, represented with Mie potentials of variable range, both with the solvent and with the other ions. A full description of the intermolecular potential requires the like energetic parameters as well as the cross-interaction parameters. The like parameters include the diameter σii, the interaction energy ϵii, and the repulsive and attractive exponents of the Mie potential, λr,ii and λa,ii, respectively. Similarly, the cross-interaction parameters include the unlike diameter σij, the unlike interaction energy ϵij, and the corresponding repulsive and attractive Mie exponents, λr,ij and λa,ij.

We also consider associating molecular ions which interact with the solvent via hydrogen bonding, and which are therefore characterised by additional association parameters: the number of sites of type a on ion i (na, i), the unlike bonding energy ϵHBab, ij, and the corresponding bonding volume Kab,ij between site a on ion i and site b on solvent j. For the development of molecular ions, we propose that the model of the molecular ion should be physically consistent with the model of the smallest neutral parent molecule giving rise to the ion. The SAFT-VR Mie model of the neutral parent molecule is used as a reference for the molecular ion and consequently the parameters of the two species, including the association parameters, will be related.

The parameterisation of the intermolecular potentials for the solvent and solute species dictate the fidelity of the proposed model. However, the complexity of our model demands the determination of a significant number of parameters, both for pure species and for unlike binary interactions. This results in a large parameter space which is known to be degenerate, yielding significant variability in the description of individual species. For charged species, the model development of electrolytes with equations of state is further complicated by the underlying premise that the ionic species can only be assessed in solution. In order to simplify the parameter estimation problem, we seek to limit the number of free parameters by assigning reasonable estimates to those parameters for which it is possible to use physical relationships to determine their value in the model. In our work, depending on the type of species i and j, the cross-interaction parameters are obtained either via combining rules or by parameter estimation, while pure-component parameters are assigned a priori values whenever possible.

3.1. Solvent model and static permittivity

We consider only aqueous electrolyte solutions and treat the solvent using the Mie model for water based on the Lennard-Jones association kernel developed in earlier work [Citation37,Citation40]. In this model, the water molecule comprises a single segment with four off-centre association sites, two of which are of the hydrogen type ‘H’ (nH, j = 2) and two of the electron lone-pair type ‘e’ (ne, j = 2), which mediate hydrogen bonding interactions. For completeness, the molecular potential parameters for this model of water are given in , together with the parameters required for the correlation of the dielectric constant (cf. Section 2.2.3).

Table 1. SAFT-VR Mie intermolecular potential parameters for the model of water (described with the effective LJ association kernel) used in our current work, taken from [Citation37], and parameters for the dielectric constant correlation, taken from [Citation21].

3.2. Ion–ion interactions

We consider first the physical geometry of the ions, addressing the diameters required to treat the ions in our model: σii and σBornii. The segment diameters of the ion models in previous work with SAFT-VRE were determined either by assigning the values of experimentally determined sizes for the bare ions [Citation13,Citation14,Citation63], or were estimated from experimental data of the thermodynamic properties of salt solutions [Citation21]. In view of the number of parameters characterising the models, we choose the former approach and assign values for σii based on experimentally reported values of the ionic sizes. With these values in hand, a standard arithmetic combining rule is used for the unlike diameter: (43)

The Born cavity diameter is commonly taken to be equal to the ionic diameter, here represented as σii, or less commonly, taken as a free parameter adjusted to provide best agreement between calculated and experimental data of the electrolyte solution. Here, we follow the work of Rashin and Honig [Citation64] who define the value of the Born cavity diameter σBornii such that the ion cavity experiences a minimum contribution from the electrons of the surrounding dielectric medium. By analysing the electron-density maps of crystals of alkali fluoride salts, Rashin and Honig proposed a 7% increase in the cavity diameter to correct for the non-sphericity of the actual ion cavity, providing physically consistent values for the Born diameters. The ion Born cavity diameters needed for the models developed in our current work are taken directly for their original paper [Citation64].

The dispersion energy ϵii between two identical ions or ϵij between any two unlike is obtained by analogy to the work of Hudson and McCoubrey [Citation21,Citation65,Citation66]. In order to obtain the dispersion energy between two ions, we relate the London [Citation67] dispersion interaction potential to the Mie intermolecular potential model given by Equation (Equation3). The London interaction potential can be expressed as a function of ionisation potentials as [Citation66] (44) where α0 is the electronic polarisability, and I is the ionisation potential of each species. In order to obtain a physical relation for the dispersive interaction energy, the London and the Mie potentials must be related. For reasons of practicality, it is easier to operate with the van der Waals integrated form of each potential ψij: (45) For the London interaction, this leads to the expression (46) while following the same procedure for the Mie potential, we obtain (47) with (48) Equating and rearranging the expressions for the integrated potentials leads to a relation for the interaction energy parameter ϵij, which can be used to estimate the value of this cross interaction for any pair of ions: (49) In this instance, the unlike attractive and repulsive Mie exponents are obtained using the combining rules derived from applying the geometric-mean criterion on the van der Waals attractive energy [Citation36]: (50) Equation (Equation50) requires knowledge of the like-ion Mie exponents, which are set a priori depending on the nature of the ion. For the atomic ions, we apply the Lennard-Jones (12-6) [Citation68] potential (which is a special case of the Mie potential with λr = 12 and λa = 6), while for the molecular ions, we adopt the form of the potential of their reference parent molecule.

3.3. Ion–solvent interactions

We consider two types of attractive interactions between the ionic and solvent species: the dispersive ion–solvent interaction which is applicable to all ions; and the hydrogen-bonding interaction between the solvent and the molecular ions possessing association sites. For the molecular ions considered here, H3O+ and OH, the smallest parent molecule is water. Consistency between these three inter-related species is achieved in part by relating their association parameters. Specifically, the ion–water association parameters are determined by scaling the association energy and bonding volume to those of the pure water–water interaction. To achieve this, the ion–water bonding volume is scaled by the corresponding unlike ion–water diameter (obtained from Equation (Equation43)), (51) and subsequently the ion–water association energy is scaled by the resulting bonding volume, (52) At this point, the ion–solvent attractive dispersion interactions remain to be determined. The exponents of the unlike Mie interaction potential are obtained using the combining rule given in Equation (Equation50), while the unlike dispersion energy between water and each ion i is treated as an adjustable parameter optimised by comparison to appropriate thermodynamic experimental data.

3.3.1. Parameter estimation

As outlined in the previous section, our proposed model development procedure requires only one adjustable parameter per ion in the case of single-solvent solutions. The estimation approach to determine the parameters makes use of experimental data for aqueous single-solute solutions, which are modelled as ternary mixtures consisting of water and the solvated ions arising from the complete dissociation of the electrolyte solute. The assumption of complete dissociation is commonly adopted in the modelling of electrolytes using EOSs; nevertheless, the physical reality of the system under consideration is known to deviate from this assumption to varying degrees [Citation69], especially at higher salt concentrations. Furthermore, as the concentration of ions in solution increases, the treatment of the solvent as a continuous dielectric medium is less appropriate. In previous work [Citation21], an upper salinity limit for the approximation of a dielectric continuum was taken as 10 molal for 1:1 electrolytes.

In order to maintain the integrity of the two aforementioned assumptions in the theory, here we choose to consider experimental data only at moderate concentrations in the estimation of the ion–water interaction energy. The consideration of experimental data for solute molalities up to 3 molal allows one to avoid biasing the ion models towards either extreme of salinity, whilst simultaneously providing a good description of the non-ideal solution behaviour at low concentrations. Aside from the careful selection of the molality range of the data sets, the temperature range considered is also restricted to a range between 278 and 473 K, so as to steer clear of the density anomaly of water close to its freezing temperature and the region close to the critical temperature of water.

The properties considered in the optimisation procedure are limited to the saturated vapour pressure p, the liquid and saturated liquid densities ρ, and the osmotic coefficient Φ of aqueous single-salt solution mixtures. We find that the use of these properties leads to robust, physically sound models for the ionic species in solution, and that other thermodynamic properties, such as the mean ionic activity coefficients of the salts, can be determined in a fully predictive manner with the resulting models. In previous work within SAFT-VRE framework [Citation21], data for the MIAC of the salts were included in the parameter estimation procedure of the ion models instead of the osmotic coefficient. Our current choice of experimental data reflects the fact that the osmotic coefficient has been studied experimentally much more extensively than the MIAC, with the latter often determined indirectly via measurements of the former using the relationship given by Equation (Equation30).

The values of the parameters are estimated by minimising an objective function consisting of the relative difference between the experimental and calculated values of the selected properties. A least-squares objective function is used following the Levenberg–Marquardt method [Citation61,Citation62]: (53) where np, o is the number of data points j for a property of type o, ωo is the weight given to property o (here we use ωo = 1 for all properties), and Xexpo, j and Xcalco, j are the experimental and calculated values of the property, respectively.

4. SAFT-VR Mie electrolyte models

The potential models of the ions developed here include nine cations and five anions, all of which can be used as constituent ions to describe multiple aqueous electrolyte solutions. We present models for five monovalent cations (Li+, Na+, K+, Rb+, H3O+), five monovalent anions (F, Cl, Br, I, and OH), and four divalent cations (Ba2 +, Ca2 +, Mg2 +, and Sr2 +).

The polarisabilities α0, i and ionisation potentials Ii required for the determination of the dispersion interactions between these ions i are readily available [Citation70–74]. For the diameters σii of the atomic ions i, we select the experimentally derived ionic diameters presented by Shannon [Citation75] corresponding to ions with a coordination number of 6 in a crystal lattice. Shannon reported a number of values of the crystal ionic diameters for a range of coordination numbers; of these, a coordination number of 6 was reported for all of the ions of interest in our current work. The choice of σii is freed from any considerations relating to the solvent environment as a direct consequence of introducing a distinct Born diameter; the latter implicitly accounts for the structure of solvent molecules around the ions. We therefore assign the size of the ions to be that of the experimentally derived crystal ionic diameter, rather than the effective ionic diameter of the solvated ion in water, as the former is expected to better represent the real size of the ions in the absence of any influence from the solvent.

In the case of the molecular hydronium cation H3O+, we refer to the SAFT-VR model of water as the closest neutral parent molecule in order to characterise its size. The protonation of water is assumed not to lead to a significant change in the size of the molecule, and hence we assign the diameter of the H3O+ ion to be that of H2O. This is considered to be a reasonable assumption given that H3O+ occurs only in aqueous solution as a solvated ion. By contrast, the hydroxyl cation OH also exists in crystalline form so, for consistency with the choice of diameters of other ions, we also assign the ionic diameter reported for OH in the work of Shannon [Citation75]. It is gratifying to find that this choice leads to a diameter value which is commensurate with the diameter of water in our SAFT-VR Mie model.

The Born cavity diameter σBornii of the atomic ions is obtained from the work of Rashin and Honig [Citation64], as discussed in Section 3.2. For the molecular ions, however, we again refer to the neutral molecule and assume that the Born diameter of both H3O+ and OH is well represented by the value for the diameter of the model of H2O, which is estimated from the bulk fluid properties of water. Rashin and Honig do not report an optimised Born diameter for H3O+, but the value they report for OH is very similar to that of our SAFT-VR Mie model of H2O, which supports our choice of approach for characterising this parameter. Helgeson and Kirkham [Citation76] have shown that there is a linear correlation between the enthalpy of solvation of ions and the inverse of the effective ionic radius, while Marcus [Citation77] has also demonstrated a correlation between the Gibbs free energy of solvation and the ionic radius. In , we apply these observations to assess our choice of Born cavity diameters for the ions and find that they correlate linearly with the experimental Gibbs free energies of solvation [Citation77–79]. This lends confidence to the values of σBornii chosen for the ions, particularly in the case of the H3O+ ion. The linear correlation of σBornii with ΔGsolv, i can be used to estimate the cavity diameter of ions as an alternative to the approach of Rashin and Honig.

Figure 2. The values of the Born cavity diameters σBornii, denoted by symbols, are shown to correlate linearly with the experimentally measured Gibbs free energy of solvation reported in [Citation77–79] for 298 K and 1.01 bar. This provides a means of validating the values assigned to σBornii in the SAFT-VR Mie models of the ions. (The dashed lines are provided as guides for the eye.)

Figure 2. The values of the Born cavity diameters σBornii, denoted by symbols, are shown to correlate linearly with the experimentally measured Gibbs free energy of solvation reported in [Citation77–79] for 298 K and 1.01 bar. This provides a means of validating the values assigned to σBornii in the SAFT-VR Mie models of the ions. (The dashed lines are provided as guides for the eye.)

The like ϵii and unlike ϵij ion–ion dispersion attractive energies are calculated using Equation (Equation49) with the polarisabilities and ionisation potentials shown in and . The ionisation potential of anions is taken to be (minus) the electron affinity of the parent species, and that of the atomic cations is taken as the higher-order ionisation potential of the parent species. These values are obtained from [Citation70] for the atomic ions, and from [Citation71–74] for the molecular ions. The resulting values for ϵii and ϵij, reported in and , respectively, follow physically reasonable trends relative to the size and charge of the ions. For atomic ions of a given charge, ϵii becomes larger with increasing ionic size. Furthermore, within a given period, the ϵii of the divalent cation is larger than that of the monovalent cation but smaller than that of the monovalent anion. The ion–ion dispersion energy is strongly dependent not only on the size of the ions but also on the form of the intermolecular potential. As well as following the correct trends, the ion–ion dispersion interactions are of a reasonable order of magnitude. This substantiates both the choice of using the Lennard-Jones potential for the atomic ions and the application of the Mie potential of the water model to represent the H3O+ and OH molecular ions.

Table 2. Values for the polarisabilities α0, i and ionisation potentials Ii of the cations i, taken from [Citation70–72].

Table 3. Values for the polarisabilities α0, i and ionisation potentials Ii of the anions i, taken from [Citation70,Citation73,Citation74].

Table 4. SAFT-VR Mie intermolecular potential parameters for the models of the solvated ions: the ion diameters σii and Born cavity diameters σBornii are obtained from the literature [Citation64,Citation75]; the like ion dispersion attractive energies ϵii are calculated using Equation (Equation49); and the ion-water unlike dispersion attractive energies are estimated using the experimental electrolyte solution data summarized in . For associating ions, the number of each site type, nH, i and ne, i, the association energy with water and the corresponding bonding volume are determined by considering the parent molecule (see section 3.3 for further details).

Table 5. Dispersion attractive energies (ϵij/kB)/K between unlike ions, calculated using Equation (Equation49).

The ion–water dispersion energy is estimated according to the procedure discussed in Section 3.3.1. The ranges of the experimental data considered are summarised in , and the sources are listed in . The estimation procedure is carried out in stages, starting first by considering all the monovalent atomic cations and anions simultaneously, using experimental solution data for 15 1:1 salts. This is followed by simultaneous estimation for all divalent atomic cations, with experimental data for 12 1:2 salts. The molecular ions are parameterised individually, using data for KOH and HBr to determine the OH–water and H3O+–water dispersion energies, respectively.

The optimal unlike parameters are shown in , and are seen to follow physically meaningful trends relative to the size and charge of the ions. The dispersion interactions between the atomic cations and water molecules increase in strength as the cations become smaller, due to the higher charge density and therefore greater polarising effect on the water molecules. The ion–water interaction of each divalent cation is also larger than that of a monovalent cation in the same period, correctly reflecting the stronger polarising effect of the smaller, higher-charge-density, divalent ions on the water molecules. For the interactions of atomic anions with water, a larger dispersion energy is obtained with increasing ionic size, as the ion becomes more polarisable.

The parameters of the molecular ions also adopt physically reasonable values, although a direct evaluation relative to the atomic ions is not possible as they differ in the range of the Mie potential and, more importantly, OH and H3O+ are modelled as associating ions. The hydronium ion is assigned three H-type sites, following the findings of MD simulations of the hydration shell of H3O+ in water by Markovitch and Agmon [Citation80]. The hydroxide ion is modelled with three e-type association sites, in line with the spectroscopic evaluation of the OH hydration shells presented by Robertson et al. [Citation81]. In aqueous solution, OH and H3O+ form hydrogen bonds with water molecules. The association parameters of these interactions are determined by scaling the hydrogen-bonding energy and bonding volume to those of pure water using Equations (Equation51) and (Equation52), based on the sizes of the molecular cores. The H3O+–water association parameters obtained with our approach are therefore the same as for pure water since the ion diameter is equal to that of water in this case.

5. Results

The adequacy of the models presented in is assessed by comparing the SAFT-VRE Mie predictions with experimental data for the saturated vapour pressure, density, osmotic coefficient, and mean ionic activity coefficient of 32 aqueous electrolyte solutions, as well as with the experimental Gibbs energy of solvation of the ions. The quality of the SAFT-VRE Mie description for these thermodynamic properties is quantified by the percentage average absolute deviation (%AAD) of each property with respect to the experimental data for that property: (54)

5.1. Description of key thermodynamic properties

Our SAFT-VRE Mie methodology provides a good description of the thermodynamic properties used in the development of the ion models within the range of thermodynamic conditions of the experimental data points used for parameter estimation. The %AAD values corresponding to the experimental dataset of are shown in . For comparison, in , we present %AAD values of the properties of the aqueous electrolyte solutions calculated with SAFT-VRE Mie from experimental data across a wide range of conditions, well beyond those considered in the model development. The expanded dataset, summarised in , includes higher salt concentrations up to 10 mol kg−1, as well as data for acid and base solutions not included in the parameter estimation procedure. The sources of these data are listed in . The high-quality performance of the SAFT-VRE Mie models is exemplified here by calculating the osmotic coefficients and densities at ambient conditions (298 K and 1.01 bar), and the vapour pressures at a range of temperatures.

Table 6. The percentage average absolute deviation %AAD of the saturated vapour pressure p, liquid density ρ, and osmotic coefficient Φ of aqueous salt solutions calculated with the SAFT-VRE Mie approach for the experimental solution data used in the parameter estimation procedure (cf. ). (The dashes indicate that experimental data for the comparison are unavailable.)

Table 7. The percentage average absolute deviation %AAD of the vapour pressure p, liquid density ρ, osmotic coefficient Φ, and mean ionic activity coefficient (MIAC) γ±, m of aqueous salt solutions calculated with the SAFT-VRE Mie approach for experimental data across a wide range of temperature and pressure conditions, subject the availability of data (cf. ). (The dashes indicate that experimental data for the comparison are unavailable.)

The description of osmotic coefficients of a range of 1:1 and 1:2 salts solutions, respectively, is illustrated in and , while the osmotic coefficients of acid and base solutions are shown in . The SAFT-VR Mie calculations can be seen to follow the trends of the experimental data, with particularly good quantitative agreement in the highly non-ideal low-molality region. Liquid-phase densities at 298 and 323 K at 1.01 bar are shown in for a selection of aqueous salt solutions. The SAFT-VR Mie representation of the density allows us to assess the methodology for the choice of diameters used to describe the ions, since this property is heavily dependent on the sizes of the species in the mixture. We specifically assess the lithium salts because Li+ is the smallest ion considered in our work, and the assumptions made regarding the ion sizes are expected to have a greater impact for the smaller ions. Given the fair agreement between the calculated densities and experimental data, we conclude that the selected crystal ionic diameters can provide a reasonable estimate of the ion size. The calculated vapour pressures of aqueous NaCl for a range of temperatures are depicted in , exemplifying the capability of the proposed model to reproduce the temperature dependence of the vapour pressure to a high level of accuracy. The results presented for these three key thermodynamic properties validate our implementation of the SAFT-VRE Mie approach for electrolytes, together with the ion models developed in our current work, as they are seen to provide a good description of the aqueous electrolyte solution properties across a broad range of conditions and compositions, with predictive capability well beyond the concentrations considered in the parameter estimation procedure.

Figure 3. The concentration dependence of the osmotic coefficient Φ for a selection of aqueous solutions of monovalent 1:1 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VRE Mie calculations, and the squares represent the experimental data obtained from the sources listed in .

Figure 3. The concentration dependence of the osmotic coefficient Φ for a selection of aqueous solutions of monovalent 1:1 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VRE Mie calculations, and the squares represent the experimental data obtained from the sources listed in Table A3.

Figure 4. The concentration dependence of the osmotic coefficient Φ for a selection of aqueous solutions of 1:2 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VRE Mie calculations, and the squares represent the experimental data obtained from the sources listed in .

Figure 4. The concentration dependence of the osmotic coefficient Φ for a selection of aqueous solutions of 1:2 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VRE Mie calculations, and the squares represent the experimental data obtained from the sources listed in Table A3.

Figure 5. The concentration dependence of the osmotic coefficient Φ for aqueous solutions of acids and bases at 298 K and 1.01 bar. The continuous curves represent the SAFT-VRE Mie calculations, and the squares represent the experimental data obtained from the sources listed in .

Figure 5. The concentration dependence of the osmotic coefficient Φ for aqueous solutions of acids and bases at 298 K and 1.01 bar. The continuous curves represent the SAFT-VRE Mie calculations, and the squares represent the experimental data obtained from the sources listed in Table A3.

Figure 6. The concentration dependence of the liquid-phase density ρ for aqueous solutions of lithium salts LiI, LiBr, and LiCl. The continuous curves and squares represent the SAFT-VRE Mie calculations and experimental data, respectively, at 298 K and 1.01 bar. The dashed curves and circles represent the SAFT-VRE Mie calculations and experimental data at 323 K and 1.01 bar. The experimental data are obtained from the sources listed in .

Figure 6. The concentration dependence of the liquid-phase density ρ for aqueous solutions of lithium salts LiI, LiBr, and LiCl. The continuous curves and squares represent the SAFT-VRE Mie calculations and experimental data, respectively, at 298 K and 1.01 bar. The dashed curves and circles represent the SAFT-VRE Mie calculations and experimental data at 323 K and 1.01 bar. The experimental data are obtained from the sources listed in Table A3.

Figure 7. The concentration dependence of the saturated vapour pressure p for aqueous solutions of NaCl for temperatures ranging from 298 to 373 K. The continuous curves represent the SAFT-VR Mie calculations, and the squares represent the experimental data obtained from the sources listed in .

Figure 7. The concentration dependence of the saturated vapour pressure p for aqueous solutions of NaCl for temperatures ranging from 298 to 373 K. The continuous curves represent the SAFT-VR Mie calculations, and the squares represent the experimental data obtained from the sources listed in Table A3.

5.2. Mean ionic activity coefficient and osmotic coefficient

A direct way of assessing the reliability of the ion models and the predictive capability of the SAFT-VRE Mie approach is via the MIAC of the aqueous salts, which is directly related to the chemical potential of the solvated ions and consequently provides a measure of how well the thermodynamic properties of the ions are represented in solution. In our current work, the MIAC is not used in the development of the ion models, and the prediction of this property can hence serve as a benchmark for ensuring that the model parameters are physically sound.

The MIAC of aqueous solution of selected salts, acids, and bases are shown in ; the %AAD of the predicted values from the corresponding experimental data for all of the salts considered are reported in . The SAFT-VRE Mie predictions for the MIAC are in excellent agreement with the experimental data despite not having been considered in the model development. By limiting the molality range of the dataset used in the parameter estimation procedure, the non-ideality of the solution at low salinity is well accounted for. As a result, this leads to very good predictions for the MIAC. Four isotherms of the MIAC of aqueous NaCl solutions are shown in for temperatures in the range 288–333 K [Citation82–88]. Correctly representing the experimental trend at low concentrations, SAFT-VRE Mie allows one to predict the decrease of the MIAC with increasing temperature for aqueous solutions of NaCl, thereby illustrating the predictive capability of the approach. At high salt concentrations, the SAFT-VRE Mie predictions continue to follow the temperature trend in the low-salinity region, which appears to be at odds with the experimental trend where the MIAC at 333 K is seen to become larger than that at 288 K; the deviations of the theoretical predictions from the experimental data are ∼20% at the highest concentrations considered. We should note, however, that there is a certain amount of scatter in the experimental data at high salinity, and that the trend for the MIAC is inconsistent with that observed for the osmotic coefficient Φ in . As the MIAC and Φ are directly related through Equation (Equation30), one would expect a similar trend with temperature, casting some doubt on the quality of the experimental data for the MIAC at higher temperatures and concentrations. On the other hand, the approximate description of the polarity of the solvent with a dielectric continuum in the SAFT-VRE approach is expected to be less adequate at very high salt concentrations.

Figure 8. The concentration dependence of the mean ionic activity coefficient γ±, m for a selection of aqueous solutions of monovalent 1:1 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in .

Figure 8. The concentration dependence of the mean ionic activity coefficient γ±, m for a selection of aqueous solutions of monovalent 1:1 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in Table A3.

Figure 9. The concentration dependence of the mean ionic activity coefficient γ±, m for a selection of aqueous solutions of 1:2 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in .

Figure 9. The concentration dependence of the mean ionic activity coefficient γ±, m for a selection of aqueous solutions of 1:2 salts at 298 K and 1.01 bar. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in Table A3.

Figure 10. The concentration dependence of the mean ionic activity coefficient γ±, m for a selection of aqueous solutions of acids and bases at 298 K and 1.01 bar. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in .

Figure 10. The concentration dependence of the mean ionic activity coefficient γ±, m for a selection of aqueous solutions of acids and bases at 298 K and 1.01 bar. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in Table A3.

Figure 11. The concentration dependence of the mean ionic activity coefficient γ±, m for aqueous solutions of NaCl at 1.01 bar for temperatures ranging from 288 to 333 K, shown both at low (top) and high (bottom) salinity. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data [Citation82–88].

Figure 11. The concentration dependence of the mean ionic activity coefficient γ±, m for aqueous solutions of NaCl at 1.01 bar for temperatures ranging from 288 to 333 K, shown both at low (top) and high (bottom) salinity. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data [Citation82–88].

Figure 12. The concentration dependence of the osmotic coefficient Φ for aqueous solutions of NaCl at 1.01 bar for temperatures ranging from 288 to 373 K. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in .

Figure 12. The concentration dependence of the osmotic coefficient Φ for aqueous solutions of NaCl at 1.01 bar for temperatures ranging from 288 to 373 K. The continuous curves represent the SAFT-VR Mie predictions, and the squares represent the experimental data obtained from the sources listed in Table A3.

5.3. Gibbs free energy of solvation

Our approach for the implementation of the Born contribution in the SAFT-VRE Mie EOS is evaluated by assessing the description of the Gibbs free energy of solvation ΔGsolv, i of the individual ions i in aqueous solution. The predictions of ΔGsolv, i are presented in , alongside the experimentally determined values obtained from [Citation77–79]. These predictions are a significant improvement over those achieved in the previous implementation of SAFT-VRE [Citation21], where the Born diameter was not differentiated from the segment diameter of the ion.

Table 8. Free energy of solvation ΔGsolv, i of ions i in aqueous solution: the SAFT-VRE Mie predictions are compared to the experimentally derived values reported in [Citation77–79].

By adhering to the appropriate definition of the Born diameter as the cavity formed by the ion in the solvent, we show that it is possible to obtain not only qualitative agreement with the trend of the solvation energies, but also good quantitative agreement with the experimental values. The level of description of the solvation effects achieved with our current implementation of SAFT-VRE Mie is similar to that of SAFT approaches in which one treats the polarity of ion–solvent interactions explicitly [Citation17,Citation28,Citation89].

5.4. Aqueous solubility of salts

In addition to the Gibbs free energy of solvation, it is also interesting to consider the limit of solubility of the salt, which can be calculated with a classical thermodynamic approach from knowledge of the activity coefficients of the ions using Equation (Equation41). This requires the activity coefficients of the ions, which are calculated using the SAFT-VRE Mie methodology, as well as the solubility product Ksp, MX of the salt MX. One way of estimating Ksp, MX is via tabulated data for the Gibbs free energies of formation ΔGf of the species using Equation (Equation42). The ΔGf of the salts and ions considered here are taken from the literature [Citation90] and are summarised in . It is important to note that these values should be used with caution as they are not direct experimental measurements and are reported as ‘best’ estimates rather than absolute quantities. Alternatively, the experimental solubility product Ksp, MXexp. can be calculated directly from experimental data for the mean ionic activity coefficient of the salt in saturated aqueous solution, by rearrangement of Equation (Equation41): (55) The solubility product of a salt obtained from Equation (Equation55) can be used with greater confidence since the data used in the calculation is specific to the salt in question. The experimental data for γexp±, m at salt saturation are obtained from [Citation91–100], and the values of Ksp calculated from Equation (Equation55) or taken from [Citation101] are presented in .

Table 9. Values used in Equation (Equation41) for the Gibbs free energies of formation of the solid salts and solvated ions obtained from [Citation90]. The of the salts correspond to the crystalline anhydrous salt at 298 K and 1.01 bar; and the of the ions correspond to the ion in aqueous solution at unit molality at 298 K and 1.01 bar.

Table 10. Values for the experimental solubility product Kspexp. at 298 K and 1.01 bar used for the calculation of the solubilities of the salts in aqueous solution.

We predict the solubility limits at 298 K and 1.01 bar for a number of aqueous salt solutions using the solubility equation with the solubility product obtained from both Equations (Equation42) and (Equation55); the results are presented in alongside the experimental solubility data [Citation70,Citation102–104]. It is immediately evident that our predicted solubilities for the most commonly studied salts are in better agreement with the reported experimental solubilities than for the less commonly encountered salts; it is possible that the tabulated reference data for the Gibbs free energy of formation and the solubility product for the more common salts are more reliable. This is supported by the fact that for the common salts, the two routes for the calculation of the solubility lead to similar predicted values. We should also point out that many of salts considered here have a solubility limit which is well above the salt concentration for which SAFT-VRE Mie is applicable. The range of application for the SAFT-VRE Mie approach is estimated to be at a maximum salt molality of about 10 mol kg−1 [Citation21], assuming a solvation shell for the ions with six coordinated water molecules. Beyond this salt concentration, the dielectric constant of the mixture can no longer be expected to be the same as that of the pure solvent, as is inherently assumed in our approach. As a consequence, it is not surprising that we find better predictions of the solubility limit for the salts with a solubility which falls within the limits of applicability of the theory. By contrast, for salts which have a solubility limit well beyond the capability of SAFT-VRE Mie, such as lithium salts, the solubility is highly overpredicted.

Table 11. The solubility limit msat for aqueous solutions of salts at 298 K and 1.01 bar. The SAFT-VRE Mie predictions are compared to the experimentally obtained values reported in [Citation70,Citation102–104]. The dashes denote that the Kspexp. for the salt is unavailable.

6. Conclusions

The SAFT-VR Mie expression for the Helmholtz free energy of a mixture is combined with the mean-spherical approximation for a non-restricted primitive model of an electrolyte solution. The proposed molecular model includes the solvent explicitly; in this instance, the solvent is water, which is modelled in the usual SAFT manner, as spherical (and non-polar) with association sites to mediate hydrogen-bonding attractive interactions in addition to repulsive and dispersive attractive interactions described with a Mie (generalised Lennard-Jones) potential. The ion–ion interactions are included via Coulombic interactions as well as repulsive and dispersive interactions (of the Mie potential form). Ion–solvent (ion–water) interactions are also incorporated; the repulsive and dispersive interactions are again taken to be of the Mie form, and the contribution to the Helmholtz free energy accounting for the charging of the ion in the solvent is included via the classical expression presented by Born. We refer to our novel generalised description of electrolyte systems as the SAFT-VRE Mie approach.

By combining literature values of ionic sizes with a well-founded physical description of the molecular interactions, the parameterisation of ionic models for the electrolyte solutions is significantly simplified. Each model requires only one ion–solvent interaction parameter, in this case the unlike dispersion attractive energy, to be adjusted using experimental data for the ionic solution. We chose to limit the range of concentration considered for the determination of this parameter to less than three molal in order to adhere to the inherent assumptions of the MSA primitive model, i.e. the representation of the solvent as a uniform dielectric medium, and the assumption of complete dissociation of the ionic species. We note the importance of a careful evaluation of the experimental data used for the determination of the interaction parameters and the range of concentrations considered. The performance of the resulting approach is shown to be accurate, even when higher molalities are considered. The thermodynamic properties of the aqueous electrolyte solutions such as the density, vapour pressure, and osmotic coefficient are reproduced well. Furthermore, the robustness and thermodynamic consistency of the approach is demonstrated by the high level of accuracy seen in the predictions of the mean ionic activity coefficient for a range of salts, in comparison to available experimental data. It is also of interest to highlight the improved agreement in the prediction of the Gibbs free energy of solvation of the ions, which is obtained by incorporating information of the cavity size required for the insertion of the charged species in the solvent. The diameter of the cavity is taken from literature sources, or correlated from a linear relationship between the solvation energy and the size of the ion when data are not available.

We demonstrate the reliability of the SAFT-VRE Mie approach in modelling the thermodynamic properties of aqueous solutions of strong electrolytes, including salts of monovalent as well as divalent ions. The models presented for the ions are shown to be robust, in the sense that the predictive capability for properties not considered in the development of the models has been confirmed; the parameters are also fully transferable to different parent salts. Importantly, the ion potential models obtained are seen to be physically realistic. This is achieved by introducing experimental information for a number of the model parameters; experimentally derived quantities such as the ionic and Born diameters are used directly, while experimental ionisation potentials, polarisabilities, and electron affinities inform ion–ion energetic parameters through the use of a theoretical relation developed in previous work [Citation66]. The number of adjustable parameters has therefore been reduced, compared with the previous formulation of SAFT-VRE SW [Citation21]. As a guiding principle, this will prove of great value in future work considering ions for which experimental data are scarce, as well as offering greater confidence in the prediction of thermodynamic properties in regions where experimental measurements are unavailable.

Data statement

Data underlying this article can be accessed on Zenodo at https://zenodo.org/record/155809, and used under the Creative Commons Attribution license.

Acknowledgments

We would like to thank T. Laffite for very useful discussions related to this work. This work was carried out as part of the Qatar Carbonates and Carbon Storage Research Centre (QCCSRC). DKE and AJH are grateful for the funding of the QCCSRC provided jointly by Qatar Petroleum, Shell, and the Qatar Science & Technology Park, for the support of a PhD studentship and research fellowship, respectively. GL is thankful to Pfizer Inc. for funding a PhD studentship. Additional support to the Molecular Systems Engineering Group from the Engineering and Physical Science Research Council (EPSRC) of the United Kingdom is also gratefully acknowledged.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The funding of this study was provided by QCCSRC; Pfizer Inc.; Qatar Petroleum; Shell; the Qatar Science and Technology Park; the Engineering and Physical Science Research Council (EPSRC) of the United Kingdom [grant number GR/T17595], [grant number GR/N35991], [grant number EP/E016340], [grant number EP/J014958]; the Joint Research Equipment Initiative (JREI) [GR/M94426]; and the Royal Society-Wolfson Foundation.

References

Appendix

Table A1. Overview of the experimental aqueous electrolyte solution data considered in the estimation procedure for the development of the SAFT-VRE Mie ion models. The ranges of temperatures T, maximum molality mmax, and the number of data points np per property (saturated vapour pressure p, density ρ, and osmotic coefficient Φ) are summarised for each salt.

Table A2. Overview of the experimental aqueous electrolyte solution data considered in the assessment of the SAFT-VRE Mie ion models. The ranges of temperatures T, maximum molality mmax, and the number of data points np per property (saturated vapour pressure p, density ρ, osmotic coefficient Φ, and mean ion activity coefficient γ±, m) are summarized for each salt.

Table A3. Sources of the experimental aqueous electrolyte solution data summarised in and .