Current Trends in Predictive Methods and Electrolyte Equations of State

Over the past 20–30 years, the development of thermodynamic models for electrolyte solutions has experienced remarkable progress. While in the first half of the 20th century, the thermodynamic models were essentially based on the continuum electrostatic approach, nowadays equations of state and predictive methods have been adapted to deal with electrolyte solutions. Given this panorama, in this mini-review, the recent advances in predictive methods and electrolyte equations of state are examined, as well as their performance in predicting activity coefficients and solid–liquid phase equilibrium data. Although this mini-review aims to shed light on the current progress in predictive methods and electrolyte equations of state, it also provides valuable references and information to the several models and theories that form the backbone of the thermodynamics of electrolyte solutions.


INTRODUCTION
Electrolyte thermodynamic models play an important role in the modeling of several natural phenomena and industrial processes. For instance, such thermodynamic models have been used to investigate chemical reactions and transport processes in natural waters, design caustic scrubbers for the removal of H 2 S(g) from gas streams, optimize purification processes in metallurgy and mining, and predict corrosion rates in pipelines, to name only a few examples. However, despite their importance in geology, biochemistry, and chemical engineering, the number of electrolyte thermodynamic models available in the literature is still quite limited in comparison to the number of thermodynamic models applicable to nonelectrolyte systems. 1−3 As a matter of fact, the development of thermodynamic models for electrolyte solutions is not so straightforward. Unlike molecular thermodynamic methods such as NRTL (nonrandom two-liquid model), 4 UNIQUAC (universal quasi-chemical model), 5 and UNIFAC (UNIQUAC functional-group activity coefficients), 6 the theoretical framework of electrolyte models must comply with additional constraints (e.g., electroneutrality) and other particularities (e.g., relative permittivity of solvents, ion−ion and ion−solvent interactions, choice of the reference state, etc.). As a consequence, some electrolyte models are only suitable for aqueous solutions under limited ranges of temperature, pressure, and concentration.
Notwithstanding the panorama described above, users of commercial simulation tools are usually confronted with a list of three to four electrolyte models that span from the extended Debye−Huckel equation to the excess Gibbs energy models (G E -models), such as ELECNRTL 7 (electrolyte NRTL model), OLI-MSE 8 (OLI-mixed solvent electrolyte), and eUNIQUAC 9 (extended UNIQUAC). Then, since the choice of the most suitable model for a simulation is not always unique, it is fundamental that the users of chemical simulators understand the fundamentals, requirements, and limitations of the electrolyte thermodynamic models, as well as the alternative approaches available.
In view of this, in this manuscript, the current trends in predictive methods and electrolyte equations of state are reviewed. In section 2, a brief discussion about the underlying concepts of G E -models and electrolyte equations of state is provided. In turn, in section 3, the use of predictive methods and equations of state in the modeling of electrolyte solutions is outlined. Although many of these models are still not available in commercial simulation tools, such models have recently grown in importance because of their predictive capabilities and potential to deal with complex mixtures over a wide range of pressure and temperature. To conclude this mini-review, the current challenges and perspectives in thermodynamic modeling of aqueous electrolyte solutions are discussed.

THEORETICAL FRAMEWORKS USED IN THE
MODELING OF ELECTROLYTE SOLUTIONS Since the early 1910s the theoretical treatment of electrolyte solutions has experienced several advances, which have resulted in the development of more realistic and comprehensive models. While in the first half of the 20th century the development of primitive models predominated, nowadays most of the electrolyte thermodynamic models rely on the estimation of excess properties and residual quantities of mixtures.
One of the key properties in chemical thermodynamics is the chemical potential μ j . Notwithstanding its conceptual relevance, the chemical potential is an abstract quantity without any equivalent in physical reality. In view of this, it is convenient to express the chemical potential in terms of auxiliary quantities that are promptly measurable, namely the fugacity and activity.
If the chemical potential is expressed in terms of the activity, one has RT x ln where μ j°i s the standard chemical potential, x j is the mole fraction of species j, R is the gas constant, T is the thermodynamic temperature, and γ j is the activity coefficient that may be estimated via the excess Gibbs energy of the electrolyte solution. Here, it is customary to assume that G E is the sum of two contributions, where G SR E and G LR E are the contributions from the short-range and long-range interactions, respectively. However, sometimes eq 2 also includes a term that accounts for middle-range interactions between ions and molecules. The long-range interactions are usually described by a Debye−Huckel type equation, while different models may be used to describe the short-range interactions, as shown in Table 1. Electrolyte models based on the exploitation of eqs 1 and 2 are called G E -models.
Alternatively, one may express the chemical potential as where φ j is the fugacity coefficient. Unlike the previous approach, here the fugacity coefficient is estimated through the relation where ΔA res is the residual Helmholtz energy, n j is the amount of substance of species j, and Z is the compressibility factor.
Electrolyte models based on the use of eqs 3 and 4 are referred to as electrolyte equations of state. 10 Overall, G E -models are best known and widely used in industrial applications. Despite that, such models require a large number of adjusting interaction parameters and their applicability is restricted to ordinary values of pressure and temperature. On the other hand, electrolyte equations of state are able to account for high values of pressure and temperature, as well as complex molecular interactions. The use of electrolyte equations of state in simulations involving mixed-solvent solutions is still under development, although concrete progress has been achieved in the last years. 11 A summary of the main electrolyte equations of state is provided in Table 2.

SYSTEMATIC REVIEW OF PREDICTIVE METHODS
AND ELECTROLYTE EQUATIONS OF STATE Next, a brief account of the ongoing advances in predictive models and electrolyte equations of state is provided. This report is not exhaustive but gives a good panorama of the most recent advances in the field. Moreover, only the main aspects of the models are considered, and further details should be obtained from the list of references.  For convenience, the performances of different models are evaluated based on the physicochemical properties and the solid−liquid equilibrium behaviors of aqueous electrolyte mixtures widely studied in the literature. To this end, theoretical data reported in some manuscripts are compared with experimental data available in refs 24 and 25. Whenever necessary, the averaging window digitalization algorithm (ΔX = 10 pixels and ΔY = 10 pixels) implemented in the software WebPlotDigitizer 26 was used to extract graphical data from 2D plots and phase diagrams.
Although this work focuses on the application of thermodynamic models to aqueous electrolyte mixtures, the predictive methods and electrolyte equations of state herein presented can also be used in the modeling of more complex systems, such as electrolyte mixtures containing amino acids 27 and alkanolamines. 28 3.1. COSMO-Based Electrolyte Models. Conductor-like screening models (COSMO) are solvation−thermodynamic methods that use the principles of quantum chemistry and statistical thermodynamics to describe molecular interactions via calculations of surface charge densities. 29 To this end, one initially considers that a molecule present in a vacuum is brought into a perfectly conducting medium of infinite relative permittivity ( Figure 1(a)). As a consequence, a polarization charge density (Figure 1 (b)) is induced at the interface between the molecule and the virtual conductor (σ-surface or simply COSMO surface), which in turn is converted into a σ-profile, a distribution function that accounts for the relative amount of molecular surface segment with charge density σ (Figure 1(c)). Next, the information provided by the σ-profile is used to estimate the σ-potential (Figure 1(d)), which describes the chemical potential of a molecular surface segment of polarity σ in the solvent s. Finally, in the last stage of the COSMO calculations, the σ-potentials of all molecular surface segments are summed to give the chemical potential of the substance in the solvent s, which is then corrected by an empirical factor related to the molecular volume and surface area in the virtual conductor environment initially created.
The theoretical framework above-discussed gave rise to two distinct predictive methods called COSMO-RS 30 (COSMO for realistic solvation) and COSMO-SAC 31 (COSMO segment activity coefficient). Then, although some refinements are still necessary, since mid-2010, COSMO-RS and COSMO-SAC have been successfully extended to electrolyte solutions.
Currently, the development of electrolyte thermodynamic models based on the COSMO-RS framework faces two main issues. 32 The first one is related to the so-called residual charge accumulated on some segments of the solvent COSMO surface. In reality, during the COSMO-RS calculations, ions are constantly interacting with the solvent molecules, so that the σ-surface of the ion might be in direct contact with the σ-surface of the solvent molecule. If this happens, the ionic σ-surface will induce an opposite polarization charge density σ A on the segments of the σ-surface of the solvent, which in turn will induce a polarization charge density σ B on the segments of the ionic σ-surface ( Figure 2). In this case, if the polarization charge density of the solvent molecule is not sufficient to completely screen the charge density of the ion, then the sum σ A + σ B is nonzero, and there will be an excess charge density on the segments of the solvent COSMO surface. As a consequence, an  energy penalty must be introduced into the calculations, which unfortunately is not taken into account in the current commercial versions of the COSMO-RS method. In addition to the issue concerning the misfit charge, the long-range interactions between ions are also not considered in the commercial software. However, such an aspect may be easily overcome, if one assumes that the ion activity coefficient is the sum of two contributions: one related to the short-range interactions, which are described by the COSMO-RS method, and another related to the long-range interactions, which may be given by the Pitzer−Debye−Huckel equation.
Given the aforementioned limitations, two approaches have been used in the development of electrolyte COSMO-RS models. In the first approach, 12 the long-range interactions are described by the Pitzer−Debye−Huckel equation, and ion− solvent clusters are considered in the COSMO-RS calculations. Here, one takes the solvation shell of an ion as an individual species and, then, assigns the respective σ-quantities through quantum chemistry calculations. By using this approach, problems might appear whenever the solution becomes very concentrated and mixed-solvents are considered. Another drawback of this approach is that the accuracy of the predictions strongly relies on the correct solvation number of every ion, but this issue was somehow solved by introducing a variable hydration number in the COSMO-RS simulations. 12 By doing so, ions are no longer treated as species with a fixed hydration number but rather as species that undergo successive (de)hydration reactions as the solution concentration changes. Through this approach, the range of validity and the accuracy of the model improved substantially, though only a few systems were tested.
Alternatively, one may also modify some terms of the equation used to estimate the ion−solvent interaction energy and treat the long-range interactions through the Pitzer− Debye−Huckel equation. Unlike the COSMO-RS-PDHS method previously described, here, the predictive capabilities of the COSMO-RS-ES model 13 might extend to highly concentrated solutions and mixed-solvents solutions, particularly if an ion pairing treatment is incorporated into the theoretical framework.
Another COSMO-based model used in the modeling of electrolyte solutions is the eCOSMO-SAC. 14 Here, the shortrange interactions are given by the COSMO-SAC model, whereas the long-range electrostatic interactions are described by the Pitzer−Debye−Huckel equation. In contrast with the COSMO-RS-PDHS and COSMO-RS-ES models, the eCOS-MO-SAC model relies on the generation of specific σ-profiles for ions or cation−anion pairs, as well as on the hypothesis of local composition. Then, according to the eCOSMO-SAC approach, the net charge rising from a distribution of ions around a solvent molecule is zero. Moreover, in very dilute solutions, ions are predominantly surrounded by solvent molecules, whereas in more concentrated solutions, ions are surrounded by solvent molecules and ions of opposite charge.
Overall, the performance of the eCOSMO-SAC model in predicting the mean ionic activity, osmotic coefficient, and vapor pressure of aqueous electrolyte solutions is good. However, the predictive capabilities of the model are severely constrained by the availability and correctness of the σ-profiles. In fact, while σprofiles for several molecules are widely available in some public and commercial repositories, only a few σ-profiles for ionic species are known in the literature. Yet, if ion-specific sigma profiles are used in the eCOSMO-SAC simulations, 14a the number of model parameters (ca. 13 segment interaction parameters plus the ionic radii) is overwhelming in comparison to other G E -models and electrolyte equations of state.

Electrolyte Group Contribution Models.
In the early 1990s, the first electrolyte thermodynamic models based on the group contribution method 33 were published. Among them, the ULPDHS (UNIFAC−Lyngby Pitzer−Debye−Huckel solvation) model proposed by Achard et al. 33b,c stands out not only by its influence on other predictive methods 12,34 but also because it is the only predictive electrolyte group contribution model implemented in a commercial process simulatorSimulis Thermodynamics.
According to the ULPDHS model, the Gibbs excess energy consists of the sum of two terms: one related to the long-range ion−ion interactions and the other related to the short-range interactions. While the long-range interactions are treated by the Pitzer−Debye−Huckel equation, the short-range interactions are described by the UNIFAC−Lingby group contribution method. In addition, Achard et al. 33b,c also incorporated the notion of average hydration number into the model in order to account for the existence of hydrated ion clusters in solution. By doing so, the description of the ion−solvent interactions is simplified and bare ions are regarded as UNIFAC independent groups.
A few years later, a new approach also based on the UNIFAC method was proposed by Yan et al. 15a However, unlike the models proposed in the early 1990s, here the excess Gibbs energy consists of the sum of three terms, where G LR E is the contribution resulting from the long-range electrostatic interactions, G SR E is the contribution from the shortrange interactions, and G MR E is the additional term due to the ion−dipole and ion−induced dipole interactions. As usual, the G LR E contribution is given by a Debye−Huckel term, whereas G SR E is estimated through the UNIFAC model. In turn, the G MR E is described by where m s is the mass of the solvent s, n k is the number of moles of the solvent group k, n i is the number of moles of the ion i, n c and n a are the number of moles of the cations c and anions a, respectively, B k,i is the interaction coefficient between the solvent group k and the ion i, B c,a is the interaction coefficient between the cation c and the anion a, and I is the ionic strength. Such interaction coefficients are related to the ionic strength through the expressions where the values of the constants b k,i , c k,i , b c,a , and c c,a are obtained through experimental data of electrolyte systems. The activity coefficients are obtained by taking the partial derivative of G E with respect to n i and n s , where MW s is the molecular weight of the solvent s and MW ms is the molecular weight of the mixed solvent that may be calculated by a proper mixing rule. Then, since eq 7 follows from a previous work authored by Li et al., 35 the thermodynamic model is referred to as LIFAC.
In the early 2000s, the original LIFAC model was revised by Mohs and Gmehling, 15b and eqs 7 and 8 were slightly modified. With such modifications, the revised LIFAC model showed more flexibility in its calculations, and the number of binary interaction parameters reduced substantially: for one main group combination, the number of parameters decreased from 11 to 6. In addition, the accuracy of the model predictions also improved, especially regarding the liquid−liquid phase equilibria. In general, both the original and modified LIFAC models provide very accurate values for the mean activity coefficients of strong electrolytes up to high concentrations ( Figure 3).

Electrolyte Equations of State.
Besides the G Emodels, the nonideality of electrolyte solutions may also be described by equations of state. In this case, the theoretical treatment of electrolyte solutions is based on the estimation of the residual Helmholtz energy ΔA res via a thermodynamic cycle similar to that shown in Figure 4, where each step shows how a specific interaction contributes to ΔA res .
Usually, the thermodynamic cycle begins with an ideal gas mixture of ions and solvent molecules, where intermolecular interactions are absent and the species are regarded as pointlike particles. Next, the ions present in the ideal gas mixture are discharged and the short-range intermolecular forces are taken into account. To this end, the Born equation and a cubic equation of state are used to estimate the contributions ΔA ⊖ Born and ΔA EoS , respectively. If specific chemical forces and structural effectshydrogen bonding interactions between the solvent molecules and/or the formation of solvation shells around the ionsare included in the model, a term ΔA Assoc needs to be incorporated in the thermodynamic cycle. Here, the Wertheim association theory is regarded as a natural choice to estimate ΔA assoc , but the statistical associating fluid theory (SAFT) may also be used. Finally, in the last two steps of the cycle the electrostatic interactions are considered: first, the ions immersed in the dielectric medium are charged, and later the long-range electrostatic forces are included in the model. For the charging process of the ions, the Born equation is used again to estimate ΔA ⊕ Born , but the permittivity of the medium must be considered. In turn, ΔA elec may be calculated by the Debye−Huckel equation or the MSA (mean spherical approximation) theory. Then, once ΔA res has been calculated via the thermodynamic cycle, one may use eq 4 to estimate the fugacity coefficient φ j and determine γ j via the relation γ j = φ j /φ j *, where φ j * is the fugacity coefficient of the species j at a reference state.
Different cubic equations of state have been used in the development of electrolyte thermodynamic models, among which stand out the Peng−Robinson equation, the Soave− Redlich−Kwong equation, and the CPA (cubic-plus-association) equation. In addition to these cubic equations of state, electrolyte models based on the SAFT framework are also numerous and well-known, as indicated in Table 2.
Basically, the CPA-based electrolyte equations of state 17,18 rely on the use of the Soave−Redlich−Kwong or the Peng− Robinson equation to calculate ΔA EoS . Moreover, the Wertheim association theory or SAFT approach is used to take into account the specific chemical interactions, whereas the longrange electrostatic interactions are described either by the MSA theory or the Debye−Huckel equation.
In turn, in the SAFT-based electrolyte equations the shortrange interactions are described through the SAFT framework, whereas the long-range electrostatic interactions are mostly given by the MSA theory, although the Debye−Huckel equation may be used as in the electrolyte PC-SAFT (perturbed-chain statistical associating fluid theory) equation proposed by Held et al. 21a Yet, in some SAFT-based models, for example, SAFT-VRE 19 (SAFT for potentials of variable range) and SAFT1-RPM 20 (heterosegmented SAFT with restricted primitive   16 Of course, changes in the thermodynamic cycle are possible and vary according to the model assumptions. model), only the self-association of water molecules is taken into account, and the influence of the solution composition on the relative permittivity of the medium is neglected. However, later Schreckenberg et al. 19c incorporated into the SAFT-VRE model a Born term, thus extending the applicability of the SAFT-VRE framework to mixed-solvent solutions too.
As alternatives to the theoretical framework described above, equations of state may also be combined with predictive methods. In this case, two approaches stand out: the VTPR-LIFAC model and the PSRK-LIFAC model. In the VTPR-LIFAC model, 22 the LIFAC model is used to provide the residual part of the molar excess Gibbs energy g E res required in the calculations of the mixture parameters a and b of the VTPR (volume-translated Peng−Robinson) equation of state, where V m is the molar volume, c is the volume translation mixture parameter, y i ′ is the salt-free mole fraction, and a ii and b ii are pure component parameters that may be obtained from critical and vapor pressure data. Hence, by coupling the LIFAC model to the VTPR equation, one may estimate salting-out effects in subcritical and supercritical fluids.
In turn, a similar procedure is used in the PSRK-LIFAC model: 23 the LIFAC method is utilized to provide the molar excess Gibbs energy g 0 E at a reference state. This quantity is essential for the calculation of the pure component parameters a ii and b i of the PSRK (predictive Soave−Redlich−Kwong) equation of state, Like the VTPR-LIFAC model, the PSRK-LIFAC model has been successfully used to estimate gas solubility in subcritical and supercritical fluids containing strong electrolytes. Although the number of works reporting the use of electrolyte equations of state has substantially increased over the past two decades, such thermodynamic frameworks are not available in the most known commercial process simulators yet. Kontogeorgis and Folas 3 provide an extensive discussion about the differences and similarities of numerous electrolyte equations of state, as well as their ranges of applicability. In general, for aqueous electrolyte mixtures under ordinary conditions of pressure and temperature, electrolyte equations of state are as successful as some of the G E -models previously presented (Figures 5 and 6). However, unlike the G E -models, electrolyte equations of state may be used to accurately estimate the solution density and treat complex mixtures at high pressures and temperatures.

SUMMARY
It is not an overstatement to say that thermodynamic models for nonelectrolyte systems are much more numerous and simple than thermodynamic models for electrolyte solutions. As a matter of fact, the development of thermodynamic models for electrolyte solutions faces several challenges that are not found in the development of models for molecular mixtures. For example, electrostatic interactions may be described either by the MSA theory or by the Debye−Huckel equation, different approaches may be used to describe the influence of the solution composition on the relative permittivity of the aqueous medium, and the interaction parameters may be ion-specific or saltspecific, to name only a few factors. As a consequence, electrolyte thermodynamics is a very diverse topic whose fundamentals are still a theme of several reviews and intense debates.
Nowadays, efforts have been mainly addressed to the development of models able to deal with extreme conditions of pressure and temperature, as well as complex mixtures of polar and associating compounds. Moreover, there is also an increasing need for models whose interaction parameters do not rely on large sets of experimental data. Concerning these aspects, predictive electrolyte models and electrolyte equations of state naturally stand out among the well-established correlative G E -models. In fact, while predictive models are particularly relevant in the chemical and pharmaceutical industries, where the thermodynamic properties of newly synthesized compounds are not always available, electrolyte equations of state have played an increasingly important role in the oil and gas industries, where high pressures and temperatures are the rule.
In addition to the points above, consortia have been created to share experimental data of electrolyte systems between academic and industrial parties and stimulate the development of reliable sets of interaction parameters for different models. A few examples of well-succeeded consortia and collaborative projects that provide thermodynamically consistent sets of interaction parameters and other experimental data for aqueous electrolyte solutions are DDBST GmbH 36a (Dortmund Data- bank), Thereda (Thermodynamic Reference Data Base), 36b KT-DTU Consortium, 36c CERE 36d (Center for Energy Resources Engineering), NIST 36e (National Institute of Standards and Technology) and its series publications, DETHERM 36f (Gesellschaft fur Chemische Technik und Biotechnologie e.V), and JESS 36g (Joint Expert Speciation System). Such organizations and projects are important partners for industries since most of the electrolyte thermodynamic models implemented in commercial process simulators are correlative. Then, as long as experimental data is available to uphold proper model parametrization, models such as ELECNRTL, OLI-MSE, and eUNIQUAC may be used in process design and other industrial tasks.
Not less important are also the initiatives toward the development and enhancement of COSMO-based electrolyte models and the implementation of new electrolyte models into computer programs. Because the availability of sigma profiles for ions and cation−anion pairs is still limited and long-range electrostatic interactions are not taken into account in the original COSMO-RS and COSMO-SAC frameworks, the application range of commercial programs based on the COSMO model is still restricted to nonelectrolyte systems. Likewise, the implementation of group predictive methods and electrolyte equations of state in process simulators is still very incipient because of the difficulties in compiling large sets of model parameters and designing cost-effective algorithms. Notwithstanding these restrictions, studies about the applicability of predictive methods and equations of state in the modeling of aqueous electrolyte solutions are more and more numerous, and soon such models will be as prevalent as the Pitzer equations and some famous G E -models.

■ ACKNOWLEDGMENTS
The author is very grateful to F. H. Quina and S. Imberti for the invitation to write this mini-review.