Saifuddin Ahmed

Nicolas Ferrando

Jean-Charles De Hemptinne
email: jean-charles.de-hemptinne@ifpen.fr

Jean-Pierre Simonin

Olivier Bernard

Olivier Baudouin

Modeling of mixed-solvent electrolyte systems

Keywords: electrolytes, mixed solvents, dielectric constant, ePPC-SAFT, salting-out

Models for mixed-solvent strong electrolytes, using an equation of state (EoS) are reviewed in this work. Through the example of ePPC-SAFT (that includes a Born term and ionic association), the meaning and the effect of each contribution to the solvation energy and the mean ionic activity coefficient are investigated. The importance of the dielectric constant is critically reviewed, with a focus on the use of a salt-concentration dependent function. The parameterization is performed using two adjustable parameters for each ion: a minimum approach distance (

) and an association energy ( ). These two parameters are optimized by fitting experimental activity coefficient and liquid density data, for all alkali halide salts simultaneously, in the range 298K to 423K. The model is subsequently tested on a large number of available experimental data, including salting out of Methane/Ethane/CO 2 /H 2 S. In all cases the deviations in bubble pressures were below 20% AADP. Predictions of vapor-liquid equilibrium of mixed solvent electrolyte systems containing methanol, ethanol are also made where deviations in bubble pressures were found to be below 10% (AADP).

Introduction

The biorefining industry involves the conversion of biomass or organic material into fuel grade biodiesel or bio-gasoline. The pre-treated biomass (feeding bio-refinery units) is a complex mixture of oxygenated hydrocarbons and water, a strongly polar solvent which forms a non-ideal mixture with the oxygenated chemicals. Water is also responsible for the degradation of processing equipment and worsening of product quality, so it needs to be separated. Aqueous solutions of salts have shown promising trend to aid in separation of these complex and oxygenated molecules encountered during the production processes of biofuel. The presence of an electrolyte causes a significant change in the equilibrium composition (especially liquid-liquid equilibrium), by altering the hydrogen bonding structure and other intermolecular forces. Hence, due to the addition of salt, the mutual solubilities change in either phase (aqueous phase and organic-rich phase). This behaviour is called the salting-out effect when the solubility decreases, and the salting-in effect when the solubility increases when adding salts [1,2]. This phenomenon is used in various industries (such as biorefining, pharmaceuticals or water treatment) for the separation of organic compounds. The use of electrolytes is however not limited to separation applications. They are often of interest in water treatment [3], geological, biological and

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

petroleum industry [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF] [START_REF] Folas | [END_REF]. Thermodynamic information is vital for various industrial processes. More specifically, chemical engineers need to know the phase equilibrium between species and enthalpies to accurately design separation processes such as distillation columns, extractive distillation, liquid-liquid extractions, etc. [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF][START_REF] Thomsen | Electrolytes Notes[END_REF][START_REF] Hendriks | [END_REF]. Due to the large number of species involved and varying operating conditions, accurate thermodynamic information is often missing, leading to inefficient design both in terms of energy efficiency and product purity. In this context, the development of a predictive thermodynamic equation of state is of primary importance [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF][START_REF] Hendriks | [END_REF][8][9][10][11][12].

Literature review

Thermodynamic models for electrolytes are numerous. The development of these models can be dated back to the 70's [13][14][15]. These models differ in applicability in terms of ranges of temperature and pressure, properties under investigation and solvents and/or salt composition. In addition to the functional form of the equation, the performance of any electrolyte model is based on many factors: the number of adjustable parameters, the accurate description of pure water (solvent), the choice of ion specific or salt specific parameters, the approach of parameterization i.e. whether salts /ion parameters are included during optimization and over what experimental data, and the study of various effects, for instance the effect of concentration on dielectric constant. While the majority of electrolyte models consider only salts that fully dissociate [16][17][18][19], a substantial minority of the models also consider partial dissociation [20][21][22][23]. A large body of work on electrolyte equations of state (EoS) has already been reviewed [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF]24,25].

Activity coefficients vs equation of state

For phase equilibrium calculations, the main property that needs to be computed is the chemical potential, which is defined as:

= , + ln (1)
where the reference state (indicated by an ) is generally taken either as the pure liquid solvent or as the fluid mixture in the ideal gas state at the same pressure and temperature. In the first case, the equation becomes

= * , + ln (2
) where i γ is the activity coefficient, which requires the use of a suitable model (e.g. Pitzer, eNRTL [26][27][28], eUNIQUAC [29]). The reference state (here indicated with *) is generally taken, for neutral molecules, at its vapor pressure ( ), while for ions, infinite dilution in the solvent, which is most often pure water. The drawback of this approach is that no pressure dependence is considered, since the activity coefficient models are generally pressure independent.

In the second case, we have: = # , , + ln (3) where i ϕ is the fugacity coefficient that requires an equation of state (EoS). The reference state (here indicated by #) is then the fluid mixture taken as an ideal gas at the same temperature T and pressure P as the fluid mixture ( : vector of composition). The logarithm of the fugacity coefficient is obtained using the mole number derivative of the volume-based residual Helmholtz energy [START_REF] Michelsen | Thermodynamic models. Fundamentals & computational aspects[END_REF]:

ln = !" # , $ !% -'% ( (4)
where ( is the compressibility factor. The relationship between the two approaches can be obtained using the definition of activity coefficients:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT = (5)
The models most often used in industry [START_REF] Lin | [END_REF] are based on activity coefficients (the most known ones are Pitzer [START_REF] Pitzer | Activity coefficients in electrolyte solutions[END_REF], modified e-NRTL [26][27][28] and e-UNIQUAC [29], as well as MSE [START_REF] Wang | [END_REF]34]). In addition to their low pressure limitation, these models rely on adjusting many interaction parameters against the available experimental data. Their predictive capability is therefore very weak. Today, many researchers try to develop electrolyte EoS [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF]35,36] which are based on the fugacity coefficient framework. They are able to account for pressure effects and allow combining the emerging statistical thermodynamic models (as the Statistical Associating Fluid Theory-SAFT-EoS based on Wertheim's association theory [37][38][39][40]) with electrolyte thermodynamics. As such, they make it possible to describe simultaneously complex molecular interactions as those occurring in bio-systems [41][42][43][44][45], for example, electrolytic effects [46][47][START_REF] Mohammad | The journal of physical chemistry[END_REF][49][50].

Yet their capabilities in modeling multi-solvent electrolyte systems are limited [START_REF] Wang | [END_REF].

Available equations of state (EoS) for mixed-solvent electrolytes

A good review on electrolyte thermodynamics was published by Loehe and Donohue [51], that included models from 1985 to 1997. Another review by Prausnitz [52] presented a brief account of electrolyte thermodynamics including their applicability in the biotechnology industry. A short review by Pinsky and Takano [START_REF] Pinsky | Computer Aided Chemical Engineering. Chapter 8: Property estimation for electrolyte systems[END_REF] presents some local composition models emphasizing computational details of activity coefficient models. Lin et al. [35] and Tan et al. [36] presented an account on electrolyte equations of state in conjunction with SAFT and electrolytic theories. In addition to these, Michelsen and Mollerup [START_REF] Michelsen | Thermodynamic models. Fundamentals & computational aspects[END_REF] presented a thorough discussion including the derivation of the Debye-Hückel theory and the theories of dipolar ions.

An extensive review of electrolyte EoSs is already presented by Kontogeorgis [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF] and Maribo-Mogensen [START_REF] Maribo-Mogensen | [END_REF], an extension of those reviews including some newer ones is presented here in table 1. This table presents the models with a mention of their parent/base equation, the type of electrolyte term employed, whether or not the Born term is used, and some indications regarding the functional form of the dielectric constant needed in the electrolyte terms (i.e. density-dependent -indicated with $; solvent concentrationdependent -indicated with %; salt concentration-dependent -indicated with )). It also states, for SAFTtype models, whether ion-ion or ion-solvent association is considered.

The table shows that the models differ on the incorporation of short-range [29,[55][56][57] and long-range forces as Debye-Hückel (DH) or mean spherical approximation(MSA) [58,59]. These models can also differ on the basis of their treatment of salt dissociation: some make use of ion-specific parameters [49], others utilize salt specific parameters [22]. Still others consider a combination of both salt and ion parameters [23,60]. Some versions of CPA [35,[START_REF] Maribo-Mogensen | [END_REF]61] and SAFT [25,[62][63][64][65][66][67][68] have been extended to model mixed-solvent electrolyte systems. However, most of them rely on fitting salt specific binary parameters which compromises with the predictive capability of the model [69].

Thermodynamic cycle: significance of each term

Thermodynamic models are often constructed using a thermodynamic cycle, where each transformation, which corresponds to a specific interaction being turned on, brings in an additive contribution to the total Gibbs energy (in the case of activity coefficient models) or Helmholtz energy (in the case of equations of state). In the case of an equation of state, the residual Helmholtz energy at given volume

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
and temperature is computed. We can then use as an example the thermodynamic cycle that is proposed by Rozmus [90]:

ΔA ,-. = ΔA /0 + ΔA 12.3 + ΔA 4..50 + ΔA 3564, + ΔA 789 + ΔA :5,; + ΔA <9=8 (6)

Figure 1 Steps in forming an electrolyte thermodynamic equation of state

The starting point is the mixture containing all species considered, but in its ideal gas state. This means that the species have no volume and no interactions. They only have kinetic energy.

Discharge

In a first step, the ions are discharged so as to yield the same mixture but without charge. The energy associated with this transformation is generally described using the Born equation [99] :

Δ" > # ? @ , $, % = - A B 4D E % ( B ?@# (7)
where is the solvation diameter of the ion, ( its charge, is the permittivity of vacuum, is the electronic charge and is the Universal Gas constant.

Repulsion and dispersion

In a second step, repulsive and attractive interactions are turned on (the species are given a volume, van der Waals, polar or hydrogen-bonding potentials). This step can be modeled with any equation of state. In fact, electrolyte equations of state have been proposed since the 90's using cubic equations of state (SRK with Fürst & Renon [19], Zuo et al. [76], Lin et al. [35], or PR with Myers et al. [22]. More recent works are based on either CPA [35] or one of the SAFT versions (see Table 1).

The Structure-forming step

Note that in this second step, the ions, though they remain present throughout the cycle, are considered in the same way as neutral molecules, which means that their only interactions are shortrange repulsion and attraction. No specific structure is created. This is obviously not the case, so the next

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
transformation that must be considered is the structure-forming of the ions: water molecules will cluster around the ions to form the so-called "solvation shell". This phenomenon is also called "hydration". Water, a polar molecule, tends to align its negative center around a cation, which forms a hydration shell.

Many molecular simulation studies exhibited this structure-forming phenomenon. As an example, figure 2 shows the radial distribution functions (rdf ) in an aqueous NaCl solution at 1 mol/kg and 298.15K between the pairs Na + /O and Cl -/O (where O is the oxygen atom of the water molecules) obtained from Monte Carlo simulations [100]. Using an integral oof the rdf, this figure also shows the number of water molecules surrounding each ion, also called the coordination numbers (CN). These coordination numbers show that in the first solvation shell there is approximately the same number of water molecules surrounding both cation (Na + ) and anion (Cl -) (6-7). The structure of water-ion mixtures is discussed at length in [101]. Fürst & Renon [19] used a specific "short range" term to describe the solvation phenomenon, but Inchekel [73] showed some inconsistencies with this term and proposed, as many others, to use the Born term instead. Yet, as we demonstrate here, the objective of the Born term is different from describing the structuring effect of ions. Several other authors [41,47,77,[START_REF] Mohammad | [END_REF]93] adjust a binary interaction parameter for the dispersion term to model this phenomenon. Several authors (e.g., Rozmus [90] and Herzog et al. [96]) use the Wertheim association term for describing this phenomenon: it allows describing the disruption of the water-water hydrogen bonds and the formation of hydration interactions in the presence of ions. This choice is also made in the present work. The short-range nature and the strong hydration interactions make the use of the SAFT association term a natural choice to model this phenomenon. The complexity comes in the selection of the number of sites on each ion. As is discussed below, we have compared two different hydration numbers: one from Bockris and Reddy [102] and other from molecular simulations to find the most suitable hydration numbers.

Electrolyte terms

The next two transformations (figure 1) are related to the presence of an electric charge on the ions: the charging effect of each single ion, and the energy related to the long-range ion-ion interactions. These two phenomena are related, and the non-primitive MSA theories [25,96] consider them explicitly. However, they are complex and rather not easy to use in process simulation context. The other existing models propose a sum of the two effects. 1. The Born term is used again to describe the energy related to the turning "on" of a point charge on each individual ion. It is identical to equation [START_REF] Hendriks | [END_REF], except for the sign and the fact that now the actual liquid solution dielectric constant is used

Δ" FGH I ? @ , $, % = A B 4D E % ( B ?@# (8)
This makes that the global Born contribution become:

" ? @ = - A B 4D J1 - 1 L E % ( B σ ?@# (9)
2. Either the Debye-Hückel (DH) [58] or Mean Spherical Approximation (MSA) [103] term is employed to describe ion-ion long-range interactions.

The long-range MSA ionic interaction theories are developed assuming ions to be point charges in a solvent described as a continuum through its dielectric constant. These theories are called "primitive models" (PM). In case of restricted primitive models the ions are treated as having identical (averaged) non-zero diameters, which allows an explicit solution of the so-called screening factor. In non-restricted PM approaches the ions are treated as having different diameters leading to an implicit model with an iterative solution [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF]. The MSA theory is based on the perturbation of polar fluids where the reference system is in the Percus-Yevick approximation using the Ornstein-Zernike equation as a specific closure.

According to some findings, it appears that MSA theory has an edge over DH theory. Galindo et al.

[77] compared the MSA and DH theories for describing short-range interactions and reported that at a higher salt concentration of NaCl, densities were more accurately represented by MSA than by DH. However, for the representation of a vapor, the performance was nearly the same. A Taylor series expansion and the comparison of the mathematical form of both these theories by Lin et al. [START_REF] Lin | Development of an equation of state for solution containing electrolytes[END_REF] showed that there are very little differences when assuming the same ion diameters. A recent comparison made by Maribo-Mogensen et al. [24] showed that the two theories gave similar results, when compared numerically in terms of screening length.

Our choice in the applications below is based on the non-restricted primitive approach. It has been recognized by many authors [START_REF] Paricaud | [END_REF]106] as a good representation of the physical phenomenon.

Some thoughts and arguments related to the choices made in this work

What is solvation?

The term solvation, also called hydration in the case of water, is used with different meanings in the literature: either it means the phenomenon related to the forming of a specific structure around an ion or that of insertion of an ion from the ideal gas to the pure solvent. From the above discussion, it is clear that the second definition, which is better defined from a theoretical point of view, is, in fact, a combination of several elementary transformations: the discharge, the cavity formation (repulsive contribution), the charging process, and the structure-forming effect around water.

The Gibbs energy of hydration from the ideal gas is rather easy to compute using the equation of state (see also Schreckenberg et al. [68]):

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT ΔN # , ?@ = N # , ?@ O -N # , ?@ I , $ = ln P ?@ O Q R S TUV RWXYZ[\ ] (10)
The values can then be compared to experimental values, as for example those given in ref. [107].

In PC-SAFT the final Helmholtz energy is the sum of various contributions such as dispersion (attraction between molecule), hard-sphere (repulsion between molecules), chain (covalent bond energy between molecules), association (hydrogen bonding energy between two associating molecules) polar (interaction energy between polar molecules), MSA (ion-ion long-range interactions) and Born (so-called solvation energy). While calculating the fugacity coefficient as in equation (10) we use the net sum of all these interactions. A look at these individual terms shows the relative contribution of each phenomenon in the computation of the Gibbs energy of hydration:

ΔN # , ?@ = ^!" # , $ !% ?@ _ GF + ^!" # , $ !% ?@ _ > #` # ?@ + ^!" # , $ !% ?@ _ FGH @ + ^!" # , $ !% ?@ _ H##?F Ha ?@ + ^!" # , $ !% ?@ _ `?bH + ^!" # , $ !% ?@ _ + ^!" # , $ !% ?@ _ ? @ -ln J # c #?bd @a L (11)
Figure 3 shows the numerical values for each of these terms at infinite dilution (the condition used for computing Gibbs energy of hydration). The value of

Q R S TUV RWXYZ[\ is provided by experiment (eg.
0.2334 at 298.15K). This term (called Z in figure 3) is clearly insignificant for the Gibbs energy of solvation. The figure shows that the contribution of the Born term is the largest in absolute value (between 450 and 200 kJ/mol). The hard chain term reflects the contribution of cavity formation and the association the structure-forming of water molecules around the ions. The other contributions are absent or negligible in infinite dilution conditions. This validates the argument that the Born contribution is essential for calculating Gibbs energy of solvation (ΔN # ) as already pointed out by [68,78]. The calculations are made using the final parameters that are used in this work (set 1 below in table 7).

Should the dielectric constant be salt-composition dependent?

Both the primitive MSA and the Born terms of the ePC-SAFT model require the dielectric constant as an input. For pure compounds, the data for the dielectric constant is abundantly available and hence several correlations exist [108][30] as a function of temperatureand sometimes density or pressure to compute their dielectric constant. In our work, we utilize the correlation proposed by Schreckenberg [68].

Yet, it has been shown experimentally that the dielectric constant is affected by salt concentration [START_REF] Maribo-Mogensen | The journal of physical chemistry[END_REF][START_REF] Shilov | The journal of physical chemistry[END_REF]. The dielectric constant of electrolyte mixtures is obtained using an experimental method known as dielectric relaxation spectroscopy. This data includes the contribution of all frequencies and electric conductivity. The dielectric constant is then calculated by extrapolating the frequency dependent permittivity to zero frequency [111]. While there is always some level of ambiguity associated to the experimental dielectric constant due to the choice of the relaxation model, the bigger problem is due to the fact that the experimental dielectric constants contain both equilibrium and non-equilibrium contribution. The non-equilibrium contribution is what is known as kinetic depolarization. Only the equilibrium contributions is required in electrolyte EoS.

Most authors of industrial electrolyte models consider that the dielectric constant to be used in the MSA and Born terms is independent of the salt concentration [START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF]. The reason for this is that the Born and MSA terms are developed in conditions of infinite dilution. They provide a corrective energy related to the charging of one ion or the ion-ion interactions. Yet, when considering the thermodynamic cycle that is presented above, it is clear that all steps bear on a system that has identical composition, temperature, and volume. Only the interactions between the species vary. Hence, when the charge is brought on the ions (charging Born term), or when the charge-charge ionic interactions are considered (MSA or DH term), the solution is not that of the infinite dilution reference state, but in terms of composition equal to that of the actual solution. This is why some authors [73,90] suggest using the salt concentration-dependent dielectric constant. Recently, Ignat and Shilov [START_REF] Shilov | The journal of physical chemistry[END_REF] used a concentration-dependent dielectric constant of several alkali halide salts to calculate their activity in an aqueous solution without any parameter adjustment which showed a semi-quantitative agreement to experimental data.

Another argument can be found in the fact that because applications generally deal with high salinities, the Born and MSA terms with infinite dilution as reference state are far beyond their application range. Hence, additional corrections would be needed anyhow. In fact, Maribo-Mogensen [START_REF] Maribo-Mogensen | The journal of physical chemistry[END_REF][START_REF] Maribo-Mogensen | The journal of physical chemistry[END_REF] has shown that the use of such a dependency provides a non-negligible improvement on the behavior of the compositional derivative of the Helmholtz energy.

An interesting method for evaluating the effect of each term on the Mean Ionic Activity Coefficient (MIAC) eq. 25-26 was provided by Inchekel et al. [73]. The combination of equations ( 3) and (4) allow writing the mean ionic activity coefficient as follows:

The ionic activity coefficient is defined as the ratio of fugacities in the mixture and in the reference state, generally pure water. ln ?@ = ln ^ ?@ ?@ * _ = ln ?@ -ln ?@ * (12)

= ^E !" l !% ?@ -E !" * l !% ?@ _ -'% J ( ( * L = E ^!" l !% ?@ - !" * l !% ?@ _ -'% J ( ( * L (13)
Where X= all contribution from the EoS. It is then possible to plot the contribution of each of these terms to the global property. Figure 4 shows such a plot, with the ePPC-SAFT model discussed in section 3 (i.e; including a composition-dependent dielectric constant).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Figure 4 Effect of the various terms on the logarithm of the mean ionic activity coefficient (MIAC) for NaCl at 298.15K. The model used is that presented in section 3 along with the parameters presented in table 6.

In the same way as already pointed out by Inchekel et al. [73], the figure shows clearly that the most significant contributions to the MIAC are the Born and the MSA terms, which need to balance out in order to obtain the well-known curved behavior with a minimum close to 1 molal (see curve "sum"). The same plot could be shown with a model not taking into account the salt dependence of the dielectric constant. In this case, no significant change is expected from the non-electrolyte terms. The behavior of the MSA term is essentially identical, but the Born contribution would necessarily be zero. This is because:

J !" !% - !" * !% L ? @ = A B 4D m ( B J 1 - 1 * L -n E % o ( o B o op ?@# q 1 B J ! !% L r,U,@ s t u (14)
Where the dielectric constant = * and its derivative with mole fraction of ion is zero. As a consequence, another contribution should balance out with the large negative value of MSA (which is necessarily so because it doesn't contain any adjustable parameter). Most often, a large dispersive energy parameter (v o ) is used such that the dispersive term is used for that purpose. Yet, an analysis of the physical significance of the above representation puts into question this choice: it can be considered that the ionic activity coefficient describes whether the ions 'likes' being in salt water better than in pure water, or not. Indeed, the fugacity, or 'escaping tendency' of the ion is proportional with the activity coefficient. When the activity coefficient is larger than 1 (positive logarithm), the ion prefers escaping from the salty water: it prefers pure water (the reference state). When, on the contrary, the activity coefficient is smaller than 1 (negative logarithm), the ions prefer to be in salty water. Now, if we consider the terms individually, we see that the MSA term, which describes the long-range electrostatic interactions between ions, always has a negative contribution. This indicates that the corresponding potential makes the ions 'like better' being surrounded by other ions: this is directly related to the decrease in screening length as discussed by Maribo-Mogensen [24]: the presence of other ions will tend to neutralize the strong local charge that is brought by the individual ion. Hence, the local electrostatic field strength will be reduced. On the other hand, the observed positive contribution of the Born term can be explained by the fact that the presence of salts reduces the dielectric constant. The energy needed to charge an ion from the vacuum is directly related to this property: the larger the dielectric constant, the lower the energy needed (the dielectric continuum makes it more comfortable for the ions). Hence, from that point of view, the ions prefer being in pure water rather than in salty water.

The other terms have a much smaller impact on the 'ionic happiness' . We see that the association (using our frame stating that association describes the forming of a hydration shell) tends to make that ions dislike salty water. This is because the water hydrogen bond network is disrupted as a consequence of their presence. The hard chain contribution is negative, because mixing spheres of different diameters creates entropy, thus increasing the stability of the mixture.

Proposed model

The SAFT equation of state

Physically, the Statistical Associating Fluid Theory (SAFT) is based on statistical perturbation theory where perturbation terms are used to correct the Helmholtz energy of a reference system.

This work propose the use of the GC-ePPC-SAFT electrolyte model by Rozmus et al. [90,113] as a basis. It is based on PC-SAFT [56] and extended to polar molecules [114][115][116] of Jog and Chapman theory for multipolar molecules [117] (" `?bH . The group contribution method [START_REF] Tamouza | Proceedings of the Fifteenth Symposium on Thermophysical Properites[END_REF] is not used in this paper. The dispersive contribution is calculated by applying perturbation to the reference fluid (which for PC-SAFT is a hard chain) based on a Lennard-Jones type pair-potential. These hard spheres are assumed to have association sites on their surface which makes it possible to calculate the Helmholtz energy contribution due to association (" H##?F ). This can also be used to calculate chain contributions by considering an infinite association strength. The electrolyte terms "

and " ? @ are also incorporated along with these terms as proposed by Inchekel et al. [73]. We will also use in this paper the correction for non-additivity of the hard sphere diameters (NAHS standing for Non-Additive Hard Sphere) that was proposed by Trinh et al. [START_REF]Prediction of phase equilibria associated with hydrotreating process of biomass by GC-PPC-SAFT Equation of State[END_REF][START_REF] Trinh | [END_REF]. The representative equation is given in eqn. 15.

A ,-. = mA /. +A 0/42; + A 12.3 + A 4..50 + A 3564, + A <9=8 + A 789 + A :5,;

The detailed mathematical expressions of the individual terms can be found in the original papers and are omitted here. Further details related to GC-PPC-SAFT can be obtained in refs [114][115][116][121][122][123][124][125][126][127]127]. Parameters for the GC-PPC-SAFT have been determined in previous papers while the parameters that are used in this work are presented here in table 3-5.

Specificities related to the GC-ePPC-SAFT model

The model used in this work is almost identical to that already discussed by Rozmus et al. [90]. Yet, since our objective was the treatment of mixed solvents, some modifications have been introduced.

The electrolyte terms

The terms that take into account the presence of charged species (MSA and Born) have not been modified [90]. Yet, the functional form of the dielectric constant is revisited. In the current work, the values of the dielectric constants for pure solvents are calculated using correlations eq. 16 developed by Schreckenberg [68]. The parameters are recalled in table 2 and for mixtures the mixing rule is given in eqn. 17.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

# = 1 + ∑ A @ $ #}#a ~xd J x U -1L (16
)

x d = ∑ % x d, #?bd @a# ∑ % #?bd @a# •%x x U = ∑ % x U, #?bd @a# ∑ % #?bd @a# (17)
In order to take into account the effect of ionic species on the dielectric constant, Pottel's [128] model is used: (eq. [18][19].

-1 = # -1 1 -€ z •• 1 + € z •• /2 ( 18
)
where # is the dielectric constant of the (mixed) solvent, " is the ionic diameter, $ is the volume of the system.

€ z •• = A d D 6 E % † " ‡ z $ ?@# (19)

The Association term

In the present approach, and as discussed in the theoretical section above, the association term is used both for describing the solvent-solvent interactions (including mixed solvents), ion-solvent and ionion short-range interactions. The association parameters are interaction parameters. Yet, most often (including in this work), they will be computed from pure component parameters using the so-called CR1 combining rule (Derawi et al. [129]).

s ˆ= †1 -‰ o ‡ ^ s + 2 _ (20)
Š s ˆ= †1 -‹ o ‡ OE Š s Š ˆ (21)
Where " ) , • Ž , Š s and Š • Ž are attributed to pure components, and ‰ o and ‹ o are adjustable interaction parameters that allows regressing the combining rules when needed. This means that regressing the ionic parameter s , the interaction is equivalent to regress the actual interaction parameter s ˆ (when ‰ o =0).

The Non-additive hard sphere diameter term

The additional term proposed by Trinh et al. [START_REF] Trinh | [END_REF] is based on the observation that the usual combining rule that is used on the diameters should be corrected with a new parameter, which is called ' o :

x o = x + x o 2 †1 -' o ‡ (22)
where the diameter d is the temperature dependent parameter as defined by Gross & Sadowski [130] :

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT x = •1 -• exp P-3 v ]• ( 23
)
where • is the sphere softness whose value is usually fixed equal to 0.12.

Parameters from previous work

The parameters utilized in the current work are mentioned in tables (3)[START_REF] Kontogeorgis | Thermodynamic models for industrial applications. From classical and advanced mixing rules to association theories[END_REF]. These parameters are taken from the previous work.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Unit and abbreviation

- Å K K - - eV D - B m σ w ε/k ε ΑΒ /v κ ΑΒ - J d > * y Q - Methanol 2

Ion parameterization procedure

The model used in this work is almost identical to that already discussed by Rozmus et al. [90] Although the results presented by Rozmus et al. [90] were of good accuracy, densities of alkali halides were not very well represented at low salinities. This is why the water parameters have been updated. It has also been pointed by several authors that the underlying model for pure solvents must be accurate for calculating their liquid densities [65,68]. It is thus important to re-parameterize the ions and see the performance of the new model for liquid densities.

The data used for the regression is presented in table 8. Mean Ionic Activity Coefficient (MIAC) eq. 24 and solution densities (c) converted to Apparent Molar Volume (AMV) using eq. 25 were used for parameterization. Their detailed description can be obtained from Rozmus et al. [90].

Mean Ionic Activity Coefficients (MIAC)

The mean ionic activity coefficient (MIAC) represented by is most commonly used to describe the electrolyte solution behaviour. It is calculated directly from the fugacity coefficient and is given as.

,~= , , O , , → 0 ™ (23) where, represents the fugacity coefficient of component, O represents the infinite dilution fugacity coefficient, ™ represents the mole fraction of water and ,~ is thus the activity coefficient in which 'm' stand for molal basis.

The Mean Ionic Activity Coefficient (MIAC) is given by:

± = † › d oe • d ž ‡ Ÿ/ d oe ›d ž (24)
where the subscripts + and -represents the cation and anion respectively. › and • are the stoichiometric coefficients (on a full dissociation basis).

Apparent Molar Volume

The apparent molar volume (AMV) of ions is defined is given as:

± = -™ ™ › + • (25)
where and ™ are the molar volumes of the solution and of pure water respectively, and › and • are the mole fractions of the cation and the anion. It is impossible to differentiate the contribution of separate ions and hence a representation of ± is used to denote the overall AMV for the solution containing ions (anions and cations). The AMV is actually the change in the volume of the solution when salts are added into it and hence can be viewed as partial molar volume for electrolyte solutions.

The parameters that can possibly be worked on, are as follows: -The ionic diameter sigma ( ): there are in fact three: one in the hard sphere term, one in the MSA term and one in the Born term. Rozmus used a single diameter for all contributions of the SAFT EoS, equal to the Pauling diameters. However, it can be argued that the diameters to be used for the MSA term [107,135,136] are not the same as the hard sphere diameter to be used for the repulsive contribution. This is why in this work we kept the Pauling diameter for the hard sphere term but regressed the MSA diameter for each ion. In order to reduce the number of parameters per ion, the diameter for the Born term was taken equal to the Pauling diameter. -The dispersive energy: while many authors use this parameter to describe the short-range interactions, we propose, as discussed above, to use the association term to describe these short-

M A N U S C R I P T A C C E P T E D

ACCEPTED MANUSCRIPT 16 range interactions. Using both dispersive and association at once would yield too many parameters so that we follow Rozmus's suggestion to consider the ionic dispersive energy to be zero.

-The Wertheim association requires two site-site interaction parameters (an energy, s , and a volume, Š s .) in addition to the choice of a number of association sites, A # . In the same way as proposed by Rozmus, we consider the association volume equal to that of water =0.044394. The association energy is regressed individually for each ion. In literature, various authors have used different A # , shown in table 6, the choice can vary widely between results from molecular dynamics calculations or experimental data (IR, spectroscopy). The objective function used is as follows

¡¢ = E £ ¤ ¥ ±,o FHbF -±,o ¦` ±,o ¦` § ¨©± opŸ B + E £ ¤ ¥ ±,o FHbF -±,o ¦` ±,o FHbF § ¨Y± opŸ B ( 26
)
where £ o is the weight of the objective function used to select the impact of the dataset on which parameters are regressed. Their values can be found in table 8. The parameters and 4.. /v for all ions were regressed simultaneously over MIAC and AMV of all alkali halide i.e. MX (where M is metal ion, and X is halide ion, for eg. NaCl).

As shown in table 6, two sets of A are compared along with the final value of the s and parameters. The values of A are smaller in set (2) as compared to set (1). In set (2) where A are taken from Bockris et al. [101], the final regressed values of s were extremely high for certain ions (see table 6). This is because the smaller the number of sites, the larger the energy must be to reach the same value for the association energy. In contrast, the final regressed values of s in set 1 are within a reasonable range. This is clearly due to a smaller number of association sites on that ion. It can be concluded that the A values from Bockris et al. [102] cannot be used.

The final parameters used in the current work are those of set 1 in table 7. The trends of the regressed diameter and association energy s are shown in fig.

(5) as compared to the Pauling diameters. It can be seen that cations and anions neatly follow a regular trend except for Li + . The reason for the high value of of Li + may be due to its small size and high charge density. The association energy s of Li + is in decreasing order which is consistent with the decreasing charge density: the higher the charge density the higher strength of association between solvent and ion. This is an expected behavior according to work of Simonin et al [140,141]. It can be noticed that the MSA diameters optimized in this work are in good agreement with the MSA diameter values recommended in this prior publication.

Regression Results

Correlation results for Alkali halide brines

The regression results of MIAC and densities are presented in fig ( 678) and in table 7 (along with data used for regression and deviations from the model). The results of MIAC show a fairly good agreement with experimental data for 20 alkali halide salts at various temperatures. The deviations are compared from the sets (1 and 2) given in table 8. The average deviations for MIAC remain at 3.95% while it is 24% for AMV from set 1 using the current model. In contrast, the deviations observed in the work of Rozmus et al. [90] were 2.9% for MIAC and 30% for AMV. The mean deviations in MIAC are 1% larger than those from Rozmus et al. [90], however that work did not take into account the AMV of fluorides which is difficult to correlate. Moreover, improvements in the current work is significant in correlating densities at low salinities and infinite dilution condition which is essential for an improved mixed solvent electrolyte model [25,68]. It can be seen that the densities of alkali halides conform fairly well at low salinities which is the capability of the model developed in our previous work [131]. The parameters for the ions that are finally used in the model are listed in table 6 (set 1). The only adjustable parameters considered here are MSA ion specific diameters and association energies. The figure 6 is grouped according to common anions to see varying trends of salt molalities and densities. For salts of chlorides, bromides and iodides, the experimental data for MIAC and calculations from the model follow an increasing trend in the order Cs + <Rb + <K + <Na + <Li + at a given concentration and at the same temperature. However, this trend is reversed in the case of fluoride salts, which is correctly captured by our model and can be seen in figure 6(d). The reason for the reversed trend of MIAC in alkali fluoride solution is the "localized hydrolysis" as explained by Robinson and Harned [142,143]: fluor acts as a base capturing protonated water thus releasing hydroxyls. Our model does not describe this chemical equilibrium phenomenon. Yet, we can try understanding the phenomenon using the decomposition of the MIAC in its contributing terms which is presented in the figure 7. This figure is drawn to see the effect of fluoride ion (for NaF), whereas fig. 4 is drawn for NaCl. The balance between the various terms is quite different in this case: the association and hard chain terms replace the Born term in positive dominance. This is the result of various factors: the s of F -is bigger as compared to that of Cl -ion resulting in a large association contribution. Hence, the trends for the Fluorides is for a large part driven by association. Physically, this is consistent with explanation of localized hydrolysis, corresponding to strong interactions between F -ion and water, which is related to the high charge density on fluoride ion. These interactions are correctly taken into account by the association term in the ePPC-SAFT model. Other authors [144] have already observed that the physical association term could be used with success to describe chemical equilibria. The calculation of AMV from the current model showed an average AAD of 30% as compared to 24% by Rozmus et al. [90]. However in the work of Rozmus et al. AMV's of fluoride salts were not included in the regression owing to their very high deviations which could have resulted in an even higher average AAD. Likewise, the computed densities that are shown in figure 8 have a very satisfactory behavior. Fig. 9 shows the effect of the temperature on the MIAC: when temperature increases, the minimum becomes shallower and the curve levels off. It shows good results for NaCl and KCl. It is also seen that the experimentally observed trend of MIAC with respect to temperature for NaCl and KCl is reproduced. The densities of the two salts NaCl and NaBr are also presented in the fig. 10, a coherent trend is captured by the model that conforms fairly well to the experimental data. It may be observed that the slope of the experimental densities flattens at higher salt molalities, whereas the model does not capture this trend. In fig. 11 multi salts (NaCl + KCl and NaCl + KBr) aqueous solution densities predictions are compared to the corresponding experimental liquid densities at 298.15K. It should be recalled here that the calculations are predictions from the model by using only ion specific parameters and no additional parameter was fitted for these systems. The model also correctly captures the variation in liquid density as a result of varying composition of two salts within the same system. The accuracy is rather good at molalities below 2 mol kg -1 . Above this salt concentration, the model steadily overestimates the densities (which is in accordance with the observations of figure 8). The global deviation remains below around 3%.

In table 8, the absolute average deviation (AAD) for MIAC and AMV are summarized for all alkali halide salts at various temperatures along with the source of experimental data.

Solvation Gibbs energies for Alkali halide brines

The Gibbs energy of solvation discussed in section 2.3.1 is an interesting property to examine. When ions are dissolved in water from the ideal gas at infinite dilution, and the same molar volume, we showed that several effects are combined resulting in a change in Gibbs energy, called the Gibbs energy of solvation. In table 9 we show the ΔN #,#Hba predicted by the current model using parameters in table 7 (set 1). The ΔN #,#Hba is defined as the sum of the ΔN # of the cation and the anion (eq. 29).

ΔN #,#Hba = ΔN #,FHa ?@ + ΔN #,H@ ?@

where ΔN #, ?@ is computed according to eqn. 11. The results are compared to result from Das et al. [25], Schreckenberg et al. [68] and Galindo et al. [77][25]. The first work used a non-primitive model (treats solvent as molecules) while the latter two consider ions in a dielectric continuum. The results presented in table 8 show that the capability of the current model to describe the solvation is qualitatively good, however, the values as those obtained by Das et al. [25] are far better. This can be due to the fact that the current model uses only two ion specific parameters while Das et al. [25] used an additional water-salt binary parameter along with two other parameters per ion. As was shown in section 2.4.1, the largest contribution to the solvation energy is the Born term. The ionic diameter used in this term has obviously a large effect on the solvation Gibbs energy. Nevertheless, it was chosen to keep this diameter equal to the Pauling diameter so as to reduce the number of adjustable parameters.

Alkanes and acid gases with brines: salting out effect in presence of organic compounds

It is well known that the solubility of organic compounds in water decreases upon addition of salts. This is called the salting-out effect. It can be presented in the form of solubility curves (figure 12) or in terms of average deviations (table 9). The systems investigated here are water + salt + methane/ethane/ CO2/ H2S. The pure component parameter for methane, ethane, CO 2, and H 2 S are listed in table 5. A binary parameter v o (for water-alkanes) and ‰ o , ‹ o (for CO 2 /H 2 S-alkanes) are regressed using the experimental data [167,168] of the solubility of methane/ethane/light gases in water using the objective function given eq. 29.

¡¢ = E ^ P FHbF -P ¦P ¦`_ ¨µS opŸ B (29)
Table 10 presents for each system studied the binary interaction parameters fitted and the average deviations obtained along with binary/ternary system parameters and experimental data.

When a salt like NaCl is added into the water + gas system, it decreases the solubility of the gas in water. This behavior is the so-called salting out phenomenon. It is difficult to capture using electrolytes EoS [66]. Our model showed a systematic underprediction of bubble pressure. It was therefore found necessary to adjust an additional parameter between the ions and the gaseous solute. Since the ions have no dispersive energy, it is impossible to use a v o for this purpose. This is why the non-additive diameter contribution developed by Trinh and coworkers [START_REF] Trinh | [END_REF] was used. This correction describes the fact that the cross-diameter is different from the arithmetic average of the two pure component diameters. It was developed in the context of hydrogen dissolution in polar compounds, using the argument that due to the polarity, the shape of the hydrogen molecule would deform. It can be argued that the presence of ions is even more disturbing than a simple dipole moment for a neutral molecule. Hence, in the presence of an ionic charge, the polarizable molecule will be distorted in such a way that the arithmetic average of the diameters may need to be modified. This can be achieved using the ' o binary interaction parameter.

By regressing the gas-ion ' o parameters, the quality of the results is improved as shown in fig 12. This parameterization is performed on the ternary vapor liquid equilibrium data of several systems given in table 10 and are given in table 10 along with average deviations, temperature range, and the experimental data used for regression. In order to have a transferable set of parameters for gas-ion pairs, the systems of a gas with all available salts are taken at once and all ion-gas parameters are regressed simultaneously. This gives gas-ion parameters that are not system specific and can be used when studying system of salts having different ions.

The values of the regressed parameters ' o between solute-ion pairs are generally small except for Li + for which the ' o was fixed to 0.5 in order to avoid computational difficulties. Such a large value would indicate that the cross diameter of the methane and Li + pair is strongly reduced compared to the arithmetic average. For the remaining ions as can be seen in the table 10 the value of ' o remains very small in absolute value. It is found to be generally positive for cations and negative for anions. The numerical values of the corrections show no clear trend, however, so that no predictive strategy could be identified. Carbon dioxide (CO 2 ) with brines has been investigated extensively due to the possibility of injection into saline aquifers principally for carbon capture and sequestration. Mineral rock containing salts dissolve in water and interact with CO 2 . In Enhanced Oil Recovery (EOR) techniques, CO2 is also viewed as a potential gas to pressurize the reservoir, where brines are naturally present. An increasing number of studies has been published, both providing new experimental data [169][170][171] and proposing thermodynamic models [90,[172][173][174] which provide accurate representations of CO 2 brine phase equilibrium. In the current work water + CO 2 + NaCl / KCl systems are investigated (table 9) with our model. It gives a fairly accurate description of the phase behavior of these systems. A cross interaction parameter ‰ o on the association energy between CO 2 and water was first fitted to obtain accurate results on salt-free system (table 10). Adding ionic species associative bonds may be created between them and CO2, according to the combining rules (equations20 and 21). Adjusting the cross-association parameters between CO2 and ions did not improve the results so it was decided to work on the ' o cross interaction parameter in the same way as for alkanes (data reference and results are shown in table 9). The salting out effect is observed as the the model fits the experimental data (figure 11). The result are shown at 313K. A change in slope in the bubble pressure is observed for 4.0 molal salt at 323.15K due to a transition from vapor-liquid to liquid-liquid region.

Another system studied in the current work is H 2 S with brine at 393.15K (figure 15). This compound also shows associative cross-interactions with ions. Yet, in the same way as for other systems a

Mixed solvent electrolytes

Since we have seen that the model is able to describe the salting-out phenomena in alkane or acid gas + brine systems, it can now be used to study vapor-liquid equilibrium of various mixed solvent salt systems. Pure water and several pure solvents such as methanol/ethanol/1-propanol have been already parameterized in a previous work [131] and parameters are recalled in table 3. Binary mixture of methanol-water were parameterized using cross interaction parameters (‰ o , ‹ o and ' o ), here three parameter were chosen to have better accuracy in the description of binary phase equilbria (without fitting ' o , deviations were close to 10% in AADP and AADY, compared to 1.8% when including this additional correction).

Several mixed solvent systems have also been studied previously by many authors [25,[START_REF] Maribo-Mogensen | [END_REF]68,[183][184][185]. However, most of them rely on fitting an additional solvent-salt/ion dispersion parameter. In this work, the results are purely predictive in nature since all parameters have been determined previously.

As seen in figure 13, the predictions from the model agree fairly well with experimental data. It is also seen that there is a systematic under-prediction of bubble temperatures at the composition of alcohols between 0.2 and 0.4. This is due to the effect of salt on the vapor-liquid equilibria of mixed solvent systems [186,187]. Figure 13-14 show the results with mixed solvent (methanol/ethanol) salt systems. The salt molality is however rather small because the salt solubility in methanol is small [188]. The blue line indicates salt-free system. Upon addition of salt the pure water boiling temperature increases, but

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

27 interestingly the mixture bubble temperature decreases. This trend can be seen in all the methanolwater-salt systems in fig. 13. The similar trend is also observed in ethanol-water-salt systems in Figure 14.

In table 10 the vapor-liquid equilibrium results of several alcohol-water-salt systems is summarized, reporting average deviations, temperature ranges and the source of experimental data. For most systems, the deviations lie, in pressure below 1% and in composition below 10%.

A quaternary system of CO2-water-methanol-NaCl is also tested predictively. The vapor-liquid equilibrium is predicted by the model by utilizing only binary parameters between methanol-water, CO2-water and CO2-ions. The average deviations in vapor pressure were 14.5% AADP, this indicated that model can predictively estimate VLE of quaternary systems with satisfactory accuracy.

Conclusion

This work addresses several issues related to the modeling of mixed solvent electrolytes. In a first part, the methodology of the construction of a SAFT-based electrolyte equation of state is discussed. More specifically, we investigate the meaning of solvation properties and propose using both ionic association and Born term to explain the physical phenomena occurring when bringing an ion from the ideal gas to the pure solvent. We also discuss the role of the dielectric constant and the effect of its functional form on the mean ionic activity coefficient. We conclude with the observation that a salinity dependent dielectric constant is a better choice for describing the individual contributing energies that come into the calculation of the MIAC.

Using the proposed analysis, the ePPC-SAFT electrolyte model initially proposed by Rozmus et al. [113] is further refined, using improved water parameters and a dielectric constant that allows working with mixed solvents. Two ion-specific parameters are fitted on activity coefficients and densities of various alkali halides.

In a final step, the model is evaluated on a larger set of VLE and density data with good success. These data include mixed salts and mixed solvent solutions. In some cases, additional binary interaction parameters had to be adjusted. For this, since the ions are considered to have negligible dispersive interactions, the algebraic average of the hard-sphere diameters was corrected using the approach proposed by Trinh et al. [START_REF] Trinh | [END_REF], with an addition adjustable parameter ' o . The alcohol-water-salt systems was subsequently modeled fairly accurately in a fully predictive way (without any parameter adjustment between alcohol and ions).

Thus, it can be concluded with reasonable certainty that the current model can be used to predict and calculate VLE of mixed-solvent electrolyte systems using solvent-water binary parameters and ionspecific parameters.

List of abbreviations. MIAC: Mean

Figure 2

 2 Figure 2 Radial distribution function (rdf ) and coordination number (CN)between Na + /O pairs and Cl -/O pairs obtained from Monte Carlo simulations [100]

Figure 3

 3 Figure 3 Contribution of Helmholtz free energy ef eg h ijk (Y-axis) for each term of PC-SAFT and for each ion at

Figure 5

 5 Figure 5 Trend of regressed parameters in comparison to Pauling diameters. (a) MSA diameters vs Pauling diameters(b) association energy ¬ f-vs Pauling diameters. ( )cations ( ) anions.

Figure 6

 6 Figure 6 MIAC of aqueous salt solutions at 298.15K for alkali halides. The lines represent calculation from model and circles represent experimental data both in the order from top to bottom ( LiX, NaX, KX, RbX, CSX) where X is the alkali halide ion except for fluorides where this order is reversed (no LiF salt not taken).

Figure 7

 7 Figure 7 Effect of the various terms on the logarithm of the mean ionic activity coefficient (MIAC) for NaF salt.

20 Figure 8

 208 Figure 8 Densities (kg/m 3 ) (Y axis) of aqueous salt solutions at 298.15K for alkali halides. The lines represent calculation from model and circles represent experimental data both in the order from bottom to top ( LiX, NaX, KX, RbX, CSX) where X is the alkali halide ion.

Figure 9

 9 Figure 9 MIAC of salt aqueous solutions at various temperatures for (a) aq. NaCl ( 273.15K 323.15K 348.15K 373.15K 423.15K) (b) aq. KCl ( 273.15K 291.15K 308.15K 323.15K) The lines represent calculation from model and circles represent experimental.

Figure 10 Figure 11

 1011 Figure 10 Liquid density(kg/m3) of salt aqueous solutions at various temperatures. The lines represent calculation from model and circles represent experimental.( 273.15K 323.15K 373.15K 423.15K 473.15K)

Figure 12

 12 Figure 12 Solubility of alkanes and acid gases in water at two different salt concentrations red(0 molal) blue (1 molal) green (4 molal) for systems (a) CH4-H2O-NaCl system at 324 K (b) C2H6-H2O-NaCl system at 353K (c) CO2-H2O-KCl system at 313K (d) H2S-H2O-NaCl system at 393K

Figure 13 Figure 14 (

 1314 Figure 13 Isobaric Vapor-liquid equilibrium of Methanol-Water-salt system at 1 bar. Lines are prediction from model and point are experimental data. (a) water-methanol-NaCl system at varying salt molalities (upto 2.56 molal) (blue line indicates salt-free systems (water-methanol) (b) water-methanol-KCl system at varying salt molalities (0.1-2 molal) in experimental data the model prediction at 2 molal concentration of salt to obtain a smooth curve at uniform salt concentration. (c) water-methanol-NaF system at varying salt molalities(0.2-0.5) in experimental data but model prediction at 0.5 molal concentration of salt to obtain a smooth curve at uniform salt concentration. (d) water-methanol-NaBr system at varying salt molalities(0.01-3.5 molal) in experimental data but model prediction at 2.0 molal concentration of salt to obtain a smooth curve at uniform salt concentration.

29 AADT:

 29 Average absolute deviation in temperature (see definition in the list of abbreviations) AADY: Average absolute deviation in vapor phase mole fraction of alcohol (see definition in the list of abbreviations)

  Ionic activity coefficient AMV: Apparent Molar Volume VLE: Vapor-Liquid Equilibrium GC-PPC-SAFT: Group Contribution Perturbed Chain Statistical Associating Fluid Theory AADP: Average absolute deviation in vapor pressure = "" absolute deviation of the alcohol in vapor phase mole fraction = Gibbs free energy of solvation. M A N U S C R I P T A C C E P T E D

31 A

 31 "ÌÍ Non-additive hard sphere term ÎÍ" Mean spherical approximation term •{ % Born term

Table 1 Review of electrolyte equations of state. This review presents an extension of the earlier reviews[4]. The electrolyte models are classified on the basis of their parent EoS model and their model for electrolytes, whether they use a term to describe solvation or not, whether they describe ion as associating species or not and on the functional form of the dielectric constant.

 1

	Model	Reference	Parent base model	Electrolyte model	use of Born term	Dielectric constant	Ion associatio n
	e-SRK	Lin et al. [35]	SRK	nRP-MSA	Born	T, V	-
	(ion specific)			and DH
		Simon et al.[70]	SRK	DH	Born	T, V, i	-
	e-(VL+DH) (ion	Zerres et al.[71]	Van Laar	DH	-	T, V [72]	-

SRK: Soave-Redlich-Kwong; PR: Peng-Robinson; CPA: Cubic Plus Association; BH: Barker Henderson, PT: Patel-Teja, SAFT-VR+DE: SAFT Variable range +Dipole and electrolyte; MSA: Mean spherical approximation, nPMSA: non-Primitive MSA, RP-MSA: restricted primitive MSA, nRP-MSA: non restricted primitive MSA; MSA(LDE): MSA low density expansion, DH: Debye Huckel T, V, i : Temperature, Volume, ion concentration; I-I: Ion-Ion , I-S : Ion-solvent

Table 2 . Correlations for calculating dielectric constants of the pure solvent.

 2

	Component Water	MW 18.02	x d 0.3777 xy z y{'	x U (K) 1403.0	%AAD 0.97
	Methanol	32.05	0.5484	1011.0	0.86
	Ethanol	46.07	0.9480	732.1	1.44

Table 3 Pure component parameters of water used in this work taken from our previous work [131]. Parameters Segment number Segment diameter Temperature dependent diameter Dispersion energy Association energy Association volume Association type

 3 This parameter (softness of diameter) is fixed as 0.12 in original PC-SAFT model, however it is made adjustable in the present version of GC-PPC-SAFT.

	Sphere Softness*	Pseudo ionization energy	Dipole Moment	Dipole fraction

*

Table 4 Pure component parameters of alkanes and alcohols used in this taken from [132] work.

 4

	Parameters	Segment number	Segment diameter	Dispersion energy	Association energy	Association volume	Association type**	Pseudo ionization energy	Dipole Moment	Dipole fraction	Quadrupole moment	Quadrupole fraction

Table 5 . Number of association sites (N s ) for ions in different electrolyte models and as calculated using molecular simulation for water-NaCl system at 298.15K

 5

			PREoS(Born	ePC-		Number of association sites A Molecular	Bockris et	This
	Ion	ePC-SAFT by Herzog et al. [96]	) by Wu and Praunitz[137]	SAFT by Lee and Kim	eSAFT-LJ by Z Liu. [83]	al.[139] Y Liu, et eEoS by	al [90] Rozmus et SAFT ePPC-	m) salt molality 1-5 (average over simulations	al.[101]	work
				[138]
	Water	2	-	4	4	4	4	-	5	4
	Li+	4	-	8	7	8	7	-	5	7
	Na+	4	6	8	7	8	7	5.64	4	7
	K+	3	-	8	7	8	7	-	3	7
	Rb+	8	-	8	8	-	7	-	3	7
	Cs+	2	-	-	10	-	7	-	4	7
	F-	0	-	-	-	-	6	-	1	6
	Cl-	0	7	7	9	6	6	6.98	1	6
	Br-	0	-	7	10	6	6	-	1	6
	I-	0	-	7	12	6	6	-	5	6

Table 6 .

 6 Final ion parameters sets, " Å are fixed equal to Pauling diameter.

		Hard
	Ion	diameter sphere		Set 1			Set 2 (Bockris et al.)[101]
	Li+	" Å 1.2	A 7	Å) 5.20577	s /v « 5748.85	A 5	Å) 0.281969	s /v « 23555.2
	Na+	1.9	7	2.74386	4962.33	5	2.4962	4611.75
	K+	2.66	7	2.61254	3869.23	4	12.398	4288.29
	Rb+	2.96	7	3.72103	2967.68	3	14.4183	3806.5
	Cs+	3.38	7	5.36451	1946.54	3	1.83048	3802.97
	F-	2.72	6	1.12055	1516.2	4	-	-
	Cl-	3.62	6	3.38508	712.074	1	1.85746	11234.5
	Br-	3.9	6	5.6969	474.425	1	-	-
	I-	4.32	6	7.70546	89.1208

Table 7 . Apparent Molar Volumes (AMV), Mean Ionic Activity Coefficients (MIAC), and densities at 298 K, experimental Data, correlations, and calculations.

 7

		Number of points and weight	Maximum molality	Deviation (%AAD) set #1		References
	Salt	AMV	MIAC
		W=0.2	W=1.0 for	AMV	MIAC	AMV	MIAC	AMV	MIAC
		for all	all
					Lithium Halides
	LiF	17	-	0.03	-	7.6	-	[147]	[148]
	LiCl	16	27	12.2	5	3.09	3.89	[149,150]	[148]
	LiBr	5	27	10	5	14.15	1.60	[151]	[148]
	LiI	4	23	3.1	3	20.39	2.7	[152]	[148]
					Sodium Halides
	NaF	24	16	0.93	1	9.85	5.85	[147]	[148]
	NaCl	12	27	6.1	5	2.84	6.21	[153,154]	[148]
	NaBr	11	27	8	5	12.59	4.05	[151,155,1 56]	[148]
	NaI	21	27	9.7	5	18.20	3.66	[157-159]	[148]
					Potassium Halides
	KF	10	27	12.7	5	57.15	2.06	[160]	[148]
	KCl	9	26	4.5	4.5	36.90	1.53	[161]	[148]
	KBr	9	27	4.4	5	40.06	2.33	[162]	[148]
	KI	11	26	8	4.5	40.89	1.93	[163]	[148]
					Rubidium Halides

Table 8 . Gibbs energy of solvation ®¯k ,k°±² for alkali halides. Experimental data compared to calculation from the model (our model vs SAFT-VR+DE[25] vs SAFT-VRE[77][25] vs SAFT-VRE[68].

 8

	Salt	Exp [107]	ePPC-SAFT (this work)	SAFT-VR+DE [25]	SAFT-VRE [77][25]	SAFT-VRE [68]
			ΔG .,2 predicted Lithium Halides	ΔG .,2 predicted	ΔG .,2 predicted	ΔG .,2 predicted
	LiF	-958	-1648.91				-1653
	LiCl	-833	-1523.42		-693	29	-1359
	LiBr	-807	-1496.19		-652	51	-1343
	LiI	-772	-1461.96		-596	65	-1299
			Sodium Halides
	NaF	-853	-1227.39		-712		-1341
	NaCl	-728	-1101.89		-653	-21	-1047
	NaBr	-702	-1074.66		-616	0.77	-1031
	NaI	-667	-1040.44		-568	11.61	-937
			Potassium Halides
	KF	-781	-1020.92		-685		-1227
	KCl	-656	-895.43		-630.	59	-933
	KBr	-630	-868.2		-593	81	-917
	KI	-595	-833.973		-548	94	-873
			Rubidium Halides
	RbF	-758	-968.612		-627	-	-1141
	RbCl	-633	-843.117		-579	-	-847

Table 10 . Summary of AAD, binary system parameters for VLE of alcohol-water-salt systems.

 10

	System	Parameter (water-alcohol)	A Q	T(K) range	Average AAD (%)	Experimental data
	(1) Methanol + (2) water + (3) LiCl	‰ o =0.0352 ‹ o =0.0025 ' o =0.0178 -	58 12	298.15 -403.15 333.15	1.83 AADP 1.71 AADY 0.76 AADT	[189] [190]
					4.54 AADY
	(3) NaF	-	25	338.95 -	0.20 AADT	[186,191]
				368.65	3.33 AADY
	(3) NaCl		11	342.05 -	0.81 AADT	[186,191]
				372.65	2.34 AADY
	(3) NaBr	-	20	338.95 -	0.28 AADT	[192]
				370.55	3.47 AADY
	(3) KCl	-	31	338.35 -	0.29 AADP	[186,191]
				372.65	7.59 AADP
	(1) Ethanol + (2) water+	OH group) ‰ o =0.11628 (water-	50	423.15 298.15 -	3.83 AADY 2.12 AADP	[193]
		‹ o =-0.9978 (water-OH group)
	(3) NaF	-	10	350.55 -	0.34 AADT	[186]
				366.15	6.07 AADY
	(3) NaBr	-	23	351.46 -	0.26 AADT	[186]
				356.45	3.19 AADY
	(3) KCl	-	12	352.2 -354	0.41 AADT	[186]
					4.22 AADY
	(3) KI	-	20	351.99 -	0.30 AADT	[186]
				356.39	4.8 AADY
	(1) CO2 + (2) Methanol + (3)Water +
		Only binary	33	313.6-395.2	14.5AADP	[194]
	(4) NaCl	parameters and
		CO2-ion (l ij )

AADP: Average absolute deviation in vapor pressure (see definition in the list of abbreviations)

Acknowledgements

This work was supported by the Tuck foundation chair on "Thermodynamics for Biofuels". Financial support from PROSIM company is also acknowledged.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

26 non-additive ' o binary interaction parameter between H 2 S and each ion was fitted on experimental saturation pressures of salted system. The current model is able to give fairly accurate description of the system as compared to experimental data.


Konstantin Akhmadeev
email: konstantin.akhmadeev@ls2n.fr

Aya Houssein

Saïd Moussaoui

Einar August Høgestøl

Ingerid Tutturen

Hanne Flinstad Harbo

Steffan Daniël Bos-Haugen

Jennifer Graves

David-Axel Laplaud

Pierre-Antoine Gourraud

Steffan D Bos-Haugen

SVM-based Tool to Detect Patients with Multiple Sclerosis Using a Commercial EMG Sensor

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Multiple sclerosis (MS) is a chronic debilitating neurological disorder that mainly affects young individuals aged between 20 and 40. As a cause of neurologic disability MS is second only to trauma, having its prevalence estimated at 2.5 million worldwide in 2014. The actual cause of MS is yet to be identified, but a complex interaction between genetic and environmental factors contributes to the risk. To date, there is no reliable method to predict MS onset or progression. Successful managing of the symptoms and attacks for MS patients highly depends on accurate and timely diagnosis as well as the possibility to measure disability progression.

Diagnostic criteria for multiple sclerosis include a number of clinical and paraclinical laboratory assessments [START_REF] Mcdonald | Recommended diagnostic criteria for multiple sclerosis: Guidelines from the International Panel on the diagnosis of multiple sclerosis[END_REF], [START_REF] Polman | Diagnostic criteria for multiple sclerosis: 2010 Revisions to the McDonald criteria[END_REF]: cerebrospinal fluid analysis, study of visual evoked potentials, electromyography analysis, neuroimaging and motor/sensory/balance function tests. The latter involves various motor tasks to be accomplished by the subject: timed 25-foot walk [START_REF] Kaufman | The significant change for the Timed 25-Foot Walk in the Multiple Sclerosis Functional Composite[END_REF], 9-hole peg test, finger-to-nose test [START_REF] Alusi | A study of tremor in multiple sclerosis[END_REF], heel-knee-shin test, finger tapping, foot tapping, etc.

The most common motor manifestations of MS are muscle fatigue, spasticity and tremor. Listed symptoms involve abnormal functioning of skeletal muscles and thus affect their activation patterns. In such cases, deviations may be revealed by analysis of limb trajectories and of involved muscles' electromyography (EMG). These measurements are proven to be efficient in different studies of MS progression [START_REF] Jonkers | Muscular effort in multiple sclerosis patients during powered wheelchair manoeuvres[END_REF]- [START_REF] Kelleher | The characterisation of gait patterns of people with multiple sclerosis[END_REF]. Thus, an EMG recording along with the inertial measurement unit (IMU) data may aid to characterise presence and severity of MS.

Common MS diagnosis approaches, as those listed above, require specific equipment, procedures and clinical expertise. A lack of them may slow down or make the diagnosis impossible, which is a common case for unpopulated areas or developing countries. A possible way to overcome these difficulties is to apply a widespread cheap acquisition system, along with unified assessment protocol and automated decision-making. As such an acquisition system we propose the MYO armband (Figure 1) developed and commercialised by Thalmic Labs Inc [START_REF]Thalmic Labs -Makers of the Myo Gesture Control Armband[END_REF]. It comprises eight EMG channels and an IMU giving acceleration, orientation and rotation speed measurements in three axes. This device is wireless, cheap, easy to use, actively supported by community and can be shipped to any location. MYO armband's default software is capable of recognising five different hand gestures, based on EMG. Also, IMU sensor provides a pointer control. In academic studies, this device was applied in sign language gesture recognition [START_REF] Abreu | Evaluating Sign Language Recognition Using the Myo Armband[END_REF] and prosthetic control [START_REF] Akhmadeev | A testing system for a real-time gesture classification using surface EMG[END_REF], [START_REF] Atasoy | 24 DOF EMG controlled hybrid actuated prosthetic hand[END_REF]. Typical signal processing pipeline in these applications consists of the following steps: windowing, feature extraction, feature projection and classification using machine learning techniques [START_REF] Phinyomark | EMG feature evaluation for improving myoelectric pattern recognition robustness[END_REF], [START_REF] Hakonen | Current state of digital signal processing in myoelectric interfaces and related applications[END_REF]. Such an approach may also be effective in an application to MS diagnostics since there is no strictly defined model of how MS affects surface EMG signals or limb trajectories. Other reasons to use machine learning techniques in this case are the dimensionality of the data and the fact that measurements are to be made with two different kinds of sensors (EMG and IMU).

In the sequel, we present the entire processing pipeline including the data acquisition, signal preprocessing, dimensionality reduction and decision-making steps. At last, experimental results involving MS patients and healthy controls are discussed.

II. MATERIALS AND METHODS

This section contains description of our equipment, test protocol, data acquisition software and data processing approaches.

A. Acquisition system

As it was already mentioned, a proposed acquisition device was an armband consisting of eight blocks connected by an elastic rubber band. Each block has a differential surface EMG electrode on the side opposed to the skin. It also contains an integrated IMU providing acceleration and orientation data. The armband hardware itself performs filtering and sampling of the signals (EMG at 200 Hz and IMU at 50Hz, 10bit quantisation) and may be connected wirelessly to a PC using Bluetooth Low Energy protocol.

A custom acquisition software with graphical user interface (figure 3) was developed to conduct the acquisition. It handled several major tasks: management of acquisition database, indication of correct armband placement and real-time signal visualisation. Back-end, or communication part of the software was developed using myo-python library [START_REF] Niklasrosenstein | Myo-python: Python bindings for the Myo SDK[END_REF], that provides a Python wrapper for original MYO libraries. Signals from all the channels were visualized during the acquisition, and then stored to a database.

B. Acquisition protocol

During tests, the MYO armband was placed either on the dominant forearm or on the lower leg, depending on the test. While on the forearm, the armband was placed so that its distal side was approximately in the middle of ulna bone and its logo-block was aligned with virtual line passing through middle fi nger (see figure 2). On the lower leg, its distal side was set in the middle of tibia, with logo-block placed laterally to it, covering tibialis anterior muscle. In function of limb size armband was loosened or tightened using special clips to assert Each subject performed seven motor function tests: timed 25 foot walk (T25FW), finger tapping (FIT), finger-to-nose (FTN, foot tapping (FT), 9 hole peg test, hell-knee-shin test and Romberg's test. Within each test, two or three trials were performed, depending on subject's performance.

C. Dataset

Only the four following tests were considered in this study: T25FW, FIT, FTN and FT. Romberg's and heel-to-knee tests were rejected because of a low signal-noise ratio. In addition, 9 hole peg test was performed on two different equipments and thus was excluded from the study. Patients and healthy control (HC) formed two classes labelled by "1" and "0" respectively. Each of the four tests was used to form a separate dataset for a separate classification problem. Within the tests, each valid trial represented a sample. Later, for each person, a classification score obtained for all trials of a selected test were averaged to obtain a final score.

D. Feature extraction

Signals from all the trials were first visually analysed and regions of interest (ROIs) were manually extracted. Next, ROI lengths were equalised by trimming longer signals. Then, from each EMG channel, twelve temporal features were extracted [START_REF] Hakonen | Current state of digital signal processing in myoelectric interfaces and related applications[END_REF]: number of zero crossings, waveform length (WL), slope sign change, Wilson amplitude (WA), root mean square, mean absolute value (MAV), integrated EMG (IEMG), signal variance and auto-regressive coefficients (4th order model).

Twelve spectral features were extracted as follows: the whole spectral band of the signal was partitioned into three sub-bands (0-20Hz, 20-50Hz, 50-100Hz), within each band mean and median frequencies (MNF and MDF) were calculated, as well as energy (EN) and dispersion of spectrum. As for IMU, only accelerometer data was taken into account, extracted features were mean and median frequencies, without band partitioning.

E. Dimensionality reduction

Overall number of sample descriptors was 24 for each of eight EMG channels plus two features for each of three accelerometer channels, giving an overall dimensionality of 198 for each sample. In order prevent overfitting during classification, PCA was applied to the feature matrix of each test. Depending on the test, from six to ten principal components were kept. Features that were most represented in first principal components are EMG features: Wilson amplitude (WA), waveform length (WL) and energy (EN) in all three sub-bands (listed in table III).

F. Classification

Support vector machine (SVM) models were learned from four datasets corresponding to each of the tests, using svmlib library [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF]. Polynomial kernel was used, best values of kernel and of regularization parameter were found using grid-search. We performed 5-fold cross validation of the result. Absolute value of SVM output is distance between the sample and separating hyperplane, while its sign represents estimated class label.

III. RESULTS AND DISCUSSION

For each data sample (single test trial), SVM model output characterises the sample score, probability of being either healthy or having MS. For each subject, scores corresponding to multiple trials of the same test were averaged. Then, sigmoidal function was applied to map the average scores onto the [0,1] interval to then trace ROC curves for each of four tests 4. Best performance obtained was for T25FW (area under ROC 95% CI [0.88, 0.99]).

The obtained scores may be interpreted differently, depending on a classification strategy. Therefore, in the situations where one single decision (for example, "no MS"/"possible MS") must be provided, the results obtained from the different tests should be combined. As an example of a possible approach, we propose a rule that consists of diagnosing a subject to have MS if 1) his most positive score is larger, by its absolute value, than his scores for the other tests; or 2) if he has at least two positive scores. Confusion matrix for that approach reveals its high sensitivity (see table IV).

IV. CONCLUSION AND FUTURE WORK

In this study we assessed the classification capabilities of MYO armband and typical EMG processing pipeline, applied to a sample of MS patients and healthy subject. Timed 25 feet walk test has shown a superior performance (in terms of Possible development of this study have to consider, first of all, the acquisition protocol. In order to compare the data acquired by two different groups, a common protocol should be established, including subject's position and trajectory of walking. Also, special care should be taken considering the intra-class variability of the EMG signals caused by different skin and body fat conditions. Necessity of EMG normalisation with respect to maximum voluntary contraction level should be investigated. No common diagnostic protocol involving MYO armband may be established until these questions are assessed. Finally, as follow-up based on this study we are investigating the possibility of predicting MS disability based on regression of SVM scores on EDSS (Expanded Disability Status Scale) [START_REF] Kurtzke | Rating neurologic impairment in multiple sclerosis: An expanded disability status scale (EDSS)[END_REF]. For that purpose, another dataset consisting only of MS patients with different levels of disability based on EDSS is currently being acquired.

Fig. 1 .

 1 Fig. 1. MYO armband by Thalmic Labs Inc.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Armband positioning on the forearm, example shown for finger tapping test.

FIT WL( 2 Fig. 4 .

 24 Fig. 4. ROC curves for four motor function tests, T25FW stands for timed 25 foot walk.