Development of a model for electrostatic contribution to the osmotic

Dec 20, 2018 - ... by the cell model is compared with that measured electrostatic interaction in our previous work [1] and it is found good match betw...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/IECR

Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Development of a Model for Electrostatic Contribution to the Osmotic Coefficient of Aqueous Electrolytes Jyoti Sahu* and Vinay A. Juvekar* Department of Chemical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai 400076, India

Ind. Eng. Chem. Res. Downloaded from pubs.acs.org by UNIV DE BARCELONA on 01/08/19. For personal use only.

S Supporting Information *

ABSTRACT: This work is the extension of our previous work [Fluid Phase Equil. 2018, 460, 57; Data in Brief 2018, 19, 485]. In this work, a predictive model has been developed for the electrostatic interaction which is based on the cell model. In this model, there are three regions around the central ion. Owing to the presence of a very high electric field (approximately 1010 V·m−1) around the central ion, water molecules are highly oriented toward the central ion and have the dielectric constant less than the dielectric constant of water in the second region (partially oriented water molecules). The third region contains the water molecules with a dielectric constant as bulk water. The electrostatic interaction, computed by the cell model, is compared with the measured electrostatic interaction in our previous work, and a good match was found between them.

1. INTRODUCTION This paper is related to our previous paper,1 in which a rationale was provided to segregate the electrostatic contribution from nonelectrostatic contribution to the osmotic coefficient for more than hundred electrolytes. The nonelectrostatic contribution to the osmotic coefficient was related to the primary hydration number of the electrolyte. The primary hydration numbers were obtained from the linear part of the osmotic coefficient−molality plot. These hydration numbers show good agreement with the experimentally measured values obtained using NMR spectroscopy. The electrostatic contribution to osmotic coefficient was computed by subtracting the nonelectrostatic contribution to osmotic coefficient from the total osmotic coefficient. It was found that the electrostatic part of the osmotic coefficient largely differs from Debye−Hückel theory. This observation provides the impetus for the development of a new model for the electrostatic contribution to the osmotic coefficient. The Debye−Hückel (D-H) theory has been used to define the electrostatic interaction since 1923. D-H theory2 has become the touchstone for the theoretical study of Coulombic system.3 According to D-H theory, the electrostatic field of the central ion is screened by the presence of an ionic cloud around it. However, the mean field nature of the D-H theory is only applicable in dilute solutions. Many models have been proposed to improve the D-H theory. For example, with the help of the concept of the potential of mean force, Kirkwood and Poirier4 have expanded the power series in terms of a charging parameter, and the coefficient of the expansion of potential of mean force was derived in a systematic way using the method of semi-invariants. Glueckauf proposed an extended version of the D-H equation for the intermediate concentration of electrolyte solutions based on cubic-root law. Pitzer5 further tested the Glueckauf equation using the virial coefficient expansion. Frank and Thompson6 and Lasaiah7 also tested the Glueckauf model by exploring © XXXX American Chemical Society

concentration ranges for which the cube-root law does not hold but transforms into an approximate square-root law. Apart from the aforementioned models, the local-composition models which have been originally developed for nonelectrolytes have been modified for electrolyte solutions by the addition of a modified D-H term to account for long-range ion−ion interactions. The local-composition models8 are based on the fact that the tendency of molecules is to show some preference in choosing their immediate neighbors in a lattice-like behavior through the presence of nonrandomness factors. The models in this category include the specific-ion interaction model by Brønsted,9,10 NRTL,11−20 UNIQUAC,21−24 UNIFAC,25 Wilson equation,26 and other hybrid models.8 Predominantly, almost all thermodynamic properties can be matched with these models so that they can be used in engineering applications under general conditions for predicting these properties in terms of temperature, pressure, and chemical composition.8 However, these models have a large number of empirical/semiempirical adjustable parameters which contain insufficient information about molecular interactions, size, and shape of particles and thermodynamic properties. As a result, when the parameters having a similar definition, obtained from one model, are used in another model, inaccurate results are obtained.27 Therefore, these models have the impression of requiring continuous upgrading.27−29 All existing electrolyte models deliver accurate results for the dilute solution. However, these models do not produce the same quality of results for higher concentration. The reason Special Issue: Vinay Juvekar Festschrift Received: Revised: Accepted: Published: A

October 9, 2018 December 20, 2018 December 20, 2018 December 20, 2018 DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

cells of cation as the central ion having radius b+ and N− cells of anion as the central ion, each having radius b−. The volume of the solution equals the sum of the volume of cells. Hence

behind this inaccuracy is that there is no direct or simple rule to extend these models for moderate concentration. Generally, these problems are solved by adding more parameters, for example, the Pitzer model. Moreover, these parameters are mostly interpreted numerically rather than physically. Varela et al. (2003)30 gave a comprehensive review of the D-H theory and its extensions, in which, they also proposed an exact theory of Coulombic systems using a combination of linearized Poisson−Boltzmann (P-B) solutions. However, few authors have tried to improve the D-H theory. Xian and Song31,32 proposed a molecular D-H theory for electrolyte solution with size symmetry in which the dielectric response of an electrolyte solution is described by a linear combination of D-H-like response-modes to calculate the mean electric potential as well as the electrostatic contribution to thermodynamic properties. According to Xian and Song, the size asymmetry of an electrolyte solution yields a harge imbalanced border zone around a solute and the dielectric response to the solute is characterized by two types of charge sources, namely, a bare solute charge and a charge distribution due to size asymmetry. These two kinds of charge sources are screened by the solvent differently. They have applied their theory for binary as well as multicomponent primitive models of electrolyte solutions. Moggia33,34 has computed the mean activity coefficients/osmotic coefficient of 3−1 electrolyte solutions at 25 °C for the maximum concentration range 1.5 to 2.2m, using the Quasi-Random Lattice (QRL) model. This model considers the fact that an ionic lattice stands for the reference electrostatic configuration which is ceaselessly disordered by molecular collisions and thermal forces occurring in the solution. According to him, the convergence of the QRL model to D-H limiting law is a notoriously difficult task for lattice models within the D-H Poisson−Boltzmann frame. In this work, we have proposed a model for predicting the contribution from electrostatic interaction to the osmotic coefficient, which is consistent with our previous work.1 In the previous work, we have abbreviated the electrostatic interaction to the osmotic coefficient as, φEe and the nonelectrostatic interaction to the osmotic coefficient as, φNE e . The focus of the present work is to develop a model which fits, φEe −m data reasonably well for electrolytes of different valence states [1:1 (e.g., NaCl, LiCl, KCl, etc.), 1:2 (e.g., Na2SO4, K2PO4, etc.)]. In the concentration range from 0 to saturation limit and temperature range from 5 to 80 °C for 1:1, 1:2, and 2:1 electrolytes. We term this model as “the cell model” since we divide the electrolyte into a set of noninteracting cells. Each cell consists of three-regions and the electrostatics of these regions is evaluated using the modified version of the Frank’s two region model.35 We confine our analysis to the region beyond the electrostatic screening limit since this region is simpler to analyze than the lower concentration region. The model parameters are estimated by fitting the model to φEe −m data for wide range of electrolytes. The attempt is made to provide the physical interpretation of the trends in the fitted parameters.

V=

4 3 4 πb+N+ = πb−3N − 3 3

(1)

Note that N+ is equal to the total number of cations in the solution minus those cations which are present as co-ions. The same is true for anion, N−. As the concentration of the electrolyte increases, the ratio of the co-ions to the central ions decreases and at sufficiently high electrolyte concentration, we can consider the entire solution to contain only central ions and counterions. In the present case, we confine ourselves to a sufficiently high electrolyte concentration, such that we can neglect co-ions. Hence N+ and N− represent the total number of cations and anions in the solution, respectively. In the present paper, we have used b = b+ = b− as a fitted parameter. Each cell is divided into three regions. The first region, surrounding the central ion is a shell containing primary hydration molecules and is devoid of ions. Water molecules in this shell are highly oriented as well as less mobile due to the very high electric field in this region. This region extends from r = 0 to r = a, where a is the minimum distance of approach between the ions. The second region of spherical shell extends from r = a to r = R* and contains the denser cloud of ions. The dielectric constant of water in this region is denoted by ϵ*. Following the Frank’s approximation, ϵ* is considered to be uniform in the region a < r < R*. The value of ϵ* is less than that of bulk water, but greater than that of water in the primary hydration layer. In the region R* < r ≤ b, the field is assumed to be low enough so that the dielectric constant of water in this region is assumed to be the same as that of the bulk water (ϵ). The model is described in brief below. The detailed derivation has been provided in the Supporting Information.

3. DEVELOPMENT OF THE MODEL FOR ELECTROSTATIC INTERACTION The electrostatic contribution to the osmotic coefficient is given by the following equation φE = −

(μw − μw0 )E M w RT ∑i mi

(2)

Here (μw−μ0w)E is the electrostatic contribution to the chemical potential of water. Using the Gibbs−Duhem equation, we can relate (μw−μ0w)E to electrostatic contribution to the chemical potential of ions of the electrolytes as dμwE = − M w ∑ mi dμi E = − M w RT ∑ mi d[ln(miγi E)] i≠w

i≠w

(3)

Here, γEi is the electrostatic contribution to the activity coefficient of ion i. Combining eq 2 and 3, φE is obtained as φE = 1 +

2. THE CELL MODEL We consider an aqueous solution containing a single electrolyte. In the cell model, the solution is divided into independent spherical cells of equal radius “b” such that each cell consists of the central ion, surrounded by counterions, co-ions and water molecules. Each cell maintains electroneutrality as a whole and therefore does not interact with the neighboring cells. The number of cells equals the number of central ions. We compute the free energy of the solution alternatively with cations and anions as the central ions. In this way, the electrolyte solution contains N+

1 ∑i mi

mi

∫0 ∑ mi d ln γi E

(4)

i≠w

γEi can be related to the free energy FEI (zie), between ion i (central ion)

of electrostatic interaction, of valency zi, and the ionic cloud surrounding it, using the following equation FIE(zie) =

RT ln γi E NA

(5)

where zie is the electric charge of the ion and NA is the Avogadro number. B

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Three stages of the system used in Frank’s model: (I) a composite system consisting of the central ion and the ionic cloud; (II) ionic cloud devoid of the central ion; (III) the central ion (nucleus) devoid of the ionic cloud. The charge distribution in the ionic cloud in state II should be identical to that in state I.

To obtain the distribution of electric field E(q,r) = −∇ϕ(q,r), first the distribution of electric potential ϕ in the region surrounding the central ion is determined by solving the linearized Poisson−Boltzmann equation.

We obtain FEI (zie) as the difference between the electrostatic energy of formation of the composite system consisting of the central ion and its cloud, FE(zie), and the sum of the electrostatic energies of formation of the individual components of the system, that is, the cloud devoid of the central ion and the central ion (nucleus) devoid of the ionic cloud. The latter two free energies are denoted as FEC(zie) and FEN(zie), respectively. Thus FIE(zie) = F E(zie) − [FCE(zie) + FNE (zie)]

ϕ 1 d ji 2 dϕ zy jjr zz = 2 r 2 dr k dr { λ* ϕ 1 d ji 2 dϕ zy jjr zz = 2 r 2 dr k dr { λ

(6)

The different states which are involved on the right-hand side of eq 6 are schematically represented in Figure 1. The columbic interaction between the central ions and the cloud is calculated by an imaginary charging process, generally either by Güntelberg charging process or Debye charging process.36 According to Onsager,37 Güntelberg and Debye charging processes give discordant results when the chemicalpotential change is calculated from ion−ion interactions using the rigorous solutions of the unlinearized P-B equation. This disagreement between the chemical-potential change calculated by Debye and Güntelberg charging processes is not due to the invalidity of either of the two charging processes; it is a symptom of the logical inconsistency intrinsic in the unlinearized P-B equation. For determining the electrostatic free energies, the Güntelberg charging process needs to be performed.36 First, FEI is written as the free energy obtained by charging the system in infinitesimal steps, starting from the uncharged state. FIE(zie) =

∫0

zie

FIE(q) dq

4πε*ε0

(11)

at r = a

and

at r = R*

dϕ dr

→0 r=b

dϕ ε* dr

(13)

=ε r = R −*

dϕ dr

r = R+*

(14)

as r → b (15)

where the boundary locations R−* and R+* belong to the inner (a ≤ r ≤ R*− ) and the outer region (R* < r ≤ b) respectively. Equation 13 is the representation of the Gauss law. The solution of eq 10 and 11 yields the distribution for the electric field for regions (a ≤ r ≤ R−*) and (R*< r ≤ b). Using the expressions of the electric fields in two situations in eq 9, and the subsequent integration would allow us to obtain FE(q). In order to obtain FC(q), the free energy of the cloud devoid of the central ion, we proceed as follows. Since the both charge distribution as well as the dielectric constants of the ionic cloud in the composite system and the system without central ion remain unchanged, by virtue of the Poisson equation, the electric fields in the two situations should be related by

(8) 1





q dϕ =− 2 dr a

ϕ|r = R −* = ϕ|r = R+*

FE(q) is obtained by integrating the free energy density 2 ε0εE(q)2 over the space occupied by the cloud. Since, the ionic cloud has two regions (a ≤ r ≤ R* and R* < r ≤ b) with different dielectric constants, the free energy integral has two parts. 1 * R* ε0ε E(q , r )2 (4πr 2) dr a 2 b 1 + ε0ε E(q , r )2 (4πr 2) dr R* 2

R* < r ≤ b

Here, cj is the molar concentration of ion j in the bulk solution. Equations 10 and 11 are coupled and can be solved together using the following boundary conditions

where, FEI (q) is the free energy of the system at electrostatic equilibrium, in which, the central ion has charge q and the accompanying ionic cloud has the net charge −q. FEI (q) can be written as

F E(q) =

(10)

where λ* and λ are the Debye lengths based on dielectric constants ϵ* and ϵ, respectively; that is, ÅÄÅ ÑÉ1/2 ÅÄÅ ÑÉ1/2 ÅÅ ε*ε RT ÑÑÑ ÅÅ ϵε RT ÑÑÑ 0 0 Å Ñ Å ÑÑ λ* = ÅÅÅ 2 and λ = ÅÅÅ 2 ÑÑ Ñ ÅÅ F ∑ z j2cj ÑÑÑ ÅÅ F ∑ z j2cj ÑÑÑ j j ÅÇÅ ÑÖÑ ÅÅÇ ÑÑÖ (12)

(7)

FIE(q) = F E(q) − [FCE(q) + FNE (q)]

a ≤ r ≤ R*

(9) C

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research 1 d 2 1 d (r E) = 2 (r 2E 0) 2 (16) r dr r dr where E is the field corresponding to the composite system and E0 is that in the cloud without the central ion. Equation 16 can be solved using the boundary conditions that q E= and E 0 = 0 at r = a 4πε*ε0a 2 (17)

Although the ionic cloud is absent, in the two regions, the dielectric constants in these two regions are different and hence the electric field is discontinuous. The following boundary conditions are applied dϕN dr

The boundary condition E0 = 0 arises due to the absence of the central ion. The solution of eq 16 is 2 q jij a zyz E (q) = E(q) − j z 4πε(r )ε0 k r {

ε* (18)

Since E(q) is already obtained in the composite solution of eq 10 and 11 along with its boundary conditions, we can obtain E0(q) by combining these equations with eq 18. Further, FC(q) is calculated using the equation identical to eq 9, except with E(q) replaced by E0(q). Since there are many common terms in E(q) and E0(q). it is simpler to calculate F(q) − FC(q) directly. ÄÅ R* ÅÅ E E F (q) − FC (q) = 2πε0ÅÅÅÅε*∫ (E2 − E 02)r 2 dr ÅÅÇ a ÉÑ ÑÑ b 2 02 2 + ε∫ (E − E )r dr ÑÑÑÑ ÑÑÖ R* (19)

1 d ijj 2 dϕN yzz z=0 jr r 2 dr jk dr z{

R* < r ≤ b

dϕN dr

=ε r = R −*

ϕN → 0

r=a

at

(22)

as

(23)

dϕN dr

r = R+*

(24)

r→∞

(25)

The solution of eq 20 and 21, is EN (a ≤ r ≤ R*) =

EN (R* < r ≤ b) =

q 4πε*ε0r 2

(26)

q 4πεε0r 2

(27)

FN(q) is obtained using eq 26 and the equation is identical to eq 9, except where E(q) is replaced by EN(q). FN(q) is given as

To obtain FN(q), the electrostatic energy of the central ion, which is devoid of the cloud, the following set of equations have been solved a ≤ r ≤ R*

q * 4πε ε0r 2

ϕN |r = R −* = ϕN |r = R+*

0

1 d jij 2 dϕN yzz z=0 jr r 2 dr jk dr z{

=−

FN (q) =

(20)

q ij 1 q ij 1 1 yz 1y jj − zz + jj − zzz b{ 8πε*ε0 k a R* { 8πεε0 k R*

(28)

Combining the expressions for F(q)−FC(q) and FN(q) in eq 8, the following equation is obtained for, FEI (q) i.e., the electrostatic free energy of interaction between the central ion with charge q and the cloud.

(21)

q q ij 1 1y jj − − zzz 4πεε0 k R* b{ 4πε*ε0R* Ä É l Å | l o ÅÅ b Ñ Ñ Ä É o o o o o o ÅÅi o o * yz ji ε zyij * yz ÅÅÅÅji ε zyij * yz ij * yzÑÑÑÑ 2(R*− b)/ λ 1 + λ ÑÑÑÑ o o o R R R R o j Å o o o j z j z j z j z Å Ñ Å Ñ o o j z j z − − + + − − − 1 1 1 1 e o j z j z j z j z Å Ñ Å Ñ o o j z j z o zz jj z j z j z ÅÅjj Ñ Å Ñ o o z j z j z b o Ñ Å Ñ * * * * o o λ { ÅÅÇk ε {k λ { k o λ { k ε {k λ {ÑÑÖ ÅÅk ÑÑ o o − 1 o o o o Å Ñ o o λ o Å Ñ o o Ç Ö o o o o o o o o o o É Ä o o o Ñ Å o o o Å b Ñ Ä É o o o Å Ñ Å Ñ o o o + 1 Å Ñ Å Ñ o o * * * * i y i y i y i y o Å Ñ Å Ñ o o R zz ji ε zyjj R zz ÅÅji ε zyjj R zz jj R zzÑÑ 2(R*− b)/ λ λ Ñ o j 2(a − R*)/ λ*Å o j Å Ñ m } j z j z o + + − + + − − + e 1 1 1 1 e j z j z j z j z Å Ñ Å Ñ j z j z o j z j z j z j z o o Å Ñ Å Ñ j z j z j z j z b o o o Å Ñ Å Ñ * * * * λ { ÅÅÇk ε {k λ { k o λ { k ε {k λ {ÑÑÖ o o Å Ñ − 1 o o o Å Ñ k o o Å Ñ λ o o o Å Ñ Öo Ç o o o o o o o o o o o o o o o o o o R* y i o o o j z o o o j z o o o zz o a − R*/ λ* R*− b / λj λ o o j o o o + 2 e e j z o o o j z o o o b z j o o o jj 1 − zz o o o o o o o λ { o o q k o n ~ − m Ä É o ÅÅ ÑÑ l b Ä É 4πε*ε0λ* o o o Å Ñ Å Ñ o o + 1 Å Ñ Å Ñ oi o o λ Ñ ÑÑ o jj1 + a yzzÅÅÅÅjijj1 − R* zyzz − ijj ε yzzjijj1 + R* zyzz + ÅÅÅÅijj ε yzzjijj1 − R* zyzz − jijj1 − R* zyzzÑÑÑÑ e 2(R*− b)/ λ oo o Ñ oj oo zz ÅÅj * zjj zz jj zzÑÑ o o b Ñ * z{ÅÅÅÅjjk * zz{ jk ε* z{jjk * o o λ λ λ λ ε λ Å Ñ o k { k o 1 − λ ÑÑÑÑ { ÅÇ k { k {ÑÖ o ÅÅÇ o oo ÑÖ o o o o o m oo ÄÅ o o o Å b Ä É o o Å Å Ñ o o 2(a − R*)/ λ*i oo o o a yzÅÅÅÅijj R* zyz ij ε yzijj R* yzz ÅÅÅÅij ε yzijj R* yzz ijj R* yzzÑÑÑÑ 2(R*− b)/ λ 1 + λ o o j o o j z j z j z −e z + Åj zj1 − z − j1 + zzÑÑ e o o j1 − * zÅÅÅjjj1 + * zzz − j * zjjj1 + o o b o o λ zz{ ÅÅÅÅÇk ε* {jjk λ zz{ jjk λ {ÅÅk λ { k ε {k λ* z{ÑÑÑÖ o o k 1− λ o o ÅÅÇ o o o on o o o o o o o o o o o o o o o o n

FIE(q) =

( (

) )

( (

( (

) )

) )

( (

D

) )

| o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o } |o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o }o o ÉÑo o o ÑÑo o o o ÑÑo o o o o ÑÑo o o o ÑÑo o o o ÑÑo o o o ÑÑo o ÑÖo o ~o o o o o o o o o o o o o o o o o o ~

(29)

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Substitution of this expression in eq 7 and subsequent integration gives FEI (zie), which upon substitution into eq 7 gives γEi as NAz 2e2 N z 2e2 ij 1 1y jj − A − zzz 8πεε0 k R* b{ 8πε*ε0R* É Ä l Å | l o ÑÑ ÅÅ b Ñ o o o ÑÉ ÅÄ o o o ÅÅi o o o o o R* yzz ji ε zyijj R* yzz ÅÅÅÅji ε zyijj R* yzz ijj R* yzzÑÑÑÑ 2(R*− b)/ λ 1 + λ ÑÑÑÑ o Åjj o o o Å o o jj zzjj1 − jj zzjj1 + o ÑÑ ÅÅjj1 − 1 e − + − − Ñ Å z j z z z o o o Ñ Å z j z z z o o zÑÑ z j z ÅÅk ε* {j z k ε* {j o Åj b Ñ o o * * o Ñ Å λ λ λ λ o o o Ñ ÅÅk 1 − Ñ Å { k { k { { k o o Ö Ç o Ñ o o λ Ñ o Å o o Ö Ç o o o o o o o o o o o o Ä É o o o Ñ Å o o o Ñ Å b É Ä o o o Ñ Å Ñ Å o o o Ñ Å 1 + Ñ Å o o o * * * * y i y i y y i i Ñ Å Ñ Å o R R R R ε ε i y i y * * * λ o o Ñ Å Ñ Å z j z j z z j j − − λ λ 2( )/ 2( )/ a R R b o j z j z Ñ Å Ñ Å z j z j z z j j o m } j z j z Ñ Å e 1 1 1 1 e + + − + + − − + Ñ Å z j z j z z j j o j z j z Ñ Å o o Ñ Å z j z j z z j j o z j z j z z j j Ñ Å o o b Ñ Å * { k ε* {k * {k * {ÑÑ o Ñ Å λ λ o o Å λ ε λ k o Ñ Å 1 − o o Å k { { k o Ö Ç Ñ Å o o o λ Ñ Å o o o Ö Ç o o o o o o o o o o o o o o o o o o R* y i o o o zz jj o o o o o o * * * j λ zz o o a − R / λ R − b / λj o o o o zz jj o o e + 2e o o o o b j 1 − zz o o o j o o o z j o o 2 2 o o o λ o NAz e o { k n ~ − m É Ä o Ñ Å l * * o Å 8πε ε0λ o b Ñ É Ä o oo Ñ Å Ñ Å Å 1 + λ ÑÑÑ Ñ Å o o o o o o jij1 + a zyzÅÅÅÅijjj1 − R* yzzz − jij ε zyzijjj1 + R* yzzz + ÅÅÅÅjij ε zyzijjj1 − R* yzzz − ijjj1 − R* yzzzÑÑÑÑe2(R*− b)/ λ ÑÑ o o o o zÅÅjj j zzÑÑ zz jj zz ÅÅj * zjj zz j * zjj o o ÑÑ b Ñ o o * * * Å λ λ Ñ Å λ λ ε ε λ o o k { k { k { Ñ Å 1 − Ñ Å { k { k { { k k o o Ö Ç Å o oo λ Ñ ÑÖ ÅÇ o o o o m o Ä o o Å oo o Å o b É Ä Å o Ñ Å o o oo o 2(a − R*)/ λ*ij a yzÅÅÅÅjij R* zzy ij ε yzjij R* zyz ÅÅÅÅij ε yzjij R* zyz jij R* zyzÑÑÑÑ 2(R*− b)/ λ 1 + λ o o o o z j j z j z Å o e 1 1 1 1 1 e + + − − + − − + − Ñ Å j z j z j z z j o zÅj j o Åj z zÑ z j zj o o b o o λ zz{ ÅÅÅÇk ε* {jjk λ zz{ jjk λ* {ÅÅÅÅjk λ* z{ k ε* {jk λ* z{ÑÑÑÖ o k o 1− λ o o o ÅÇ o o n o o o o o o o o o o o o o o o o n

RT ln γi E =

( (

) )

( (

( (

) )

) )

( (

) )

| o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o } o | o o o o o o o o o o o o o o o o o o o o oo o o o o oo o }o o o ÉÑo o o ÑÑo o o ÑÑo o o oo o ÑÑo o ÑÑo o o oo o ÑÑo o o ÑÑo o o ÑÖo o ~o o o o o o o o o o o o o o o o o ~

(30)

Substituting the expression of ln γEi from eq 30 into eq 6 and integrating the expression from 0 to m and noting that φEe = 0 at m = 0, we get the following expression for φEe φeE = −

ij 1 jjIf (I ) − 4π ∑i ≠ w mi k

∫0

I

y f (I )dI zzz {

(31)

where

É ÄÅ | l ÑÑ ÅÅ b Ñ o o ÑÉÑ ÅÄÅ o o 1 + Ñ ÅÅi o o y i y i i y y * * * * Ñ Å o o i y i y R R R R ε ε *− b)/ λ λ Ñ j z j z z j j z Ñ Ñ Å Å R 2( o o z j j z z j z j j z j z Ñ Ñ Å Å o oÅÅjj1 − z j j z 1 1 1 e − + + − − − z j z j j z z Ñ Ñ Å o o j z j z j z j z j z z Ñ Ñ Å o o b Ñ Ñ Å Å * * * * o o λ { ÅÅÇk ε {k λ { k λ { k ε {k λ {ÑÑÖ ÑÑ o oÅÅÅk 1 − o o Ñ o o λ Ñ Å o o Ö Ç o o o o o Ä Éo o o Å Ñ o o Å É Ä b Ñ o o Å Ñ Ñ Å o o 1 + Å Ñ Ñ Å o o i y i i y i y y * * * * Å Ñ Ñ Å o o 2(a − R*)/ λ*ÅÅjj R zz ji ε zyjj R zz ÅÅji ε zyjj R zz jj R zzÑÑ 2(R*− b)/ λ λ Ñ Ñ } m j z j z e 1 1 1 1 e + + − + + − − + j z j z j z j z Å Ñ Ñ Å j * zj j * zj o o j z z j z z Å Ñ Ñ Å b o o Å Ñ Ñ Å * * λ { ÅÅÇk ε {k λ { k λ { k ε {k λ {ÑÑÖ o o Å Ñ 1 − o o k Å Ñ o o Å Ñ λ o o Å Ñ Ç Ö o o o o o o o o o o * y o o R i o o z j o o z j o o */ λ* R*− b / λj z − λ a R o o zz jj o o+ 2e e o o z j o o b z j o o z j 1 − o o z j o o ij y λ Fe { k zz ~ n j zz f (I ) = jjj Ä É z Å Ñ l * * ÅÅ É Ä b Ñ Ñ Ñ Å o k λ ε ε0RT { o 1 + λ ÑÑ Å Ñ Å o o ÑÑ o jij1 + a zyzÅÅÅÅijjj1 − R* yzzz − jij ε zyzijjj1 + R* zzzy + ÅÅÅÅjij ε zyzijjj1 − R* yzzz − ijjj1 − R* yzzzÑÑÑÑe2(R*− b)/ λ o Ñ o j zÅÅj j * zj j * zj z z j z z Ñ Å o b Ñ Ñ ÑÑ Å * * * o λ { ÅÇÅk ε {k λ { k λ {ÅÅk λ { k ε {k λ {ÑÑÖ k o 1 − ÑÑ o Å o λ Å ÑÖ o Ç m Ä o Å o b o ÅÄ ÅÅ ÑÉ o o o i a zyÅÅÅÅijj R* yzz ji ε zyijj R* yzz ÅÅÅÅji ε zyijj R* yzz ijj R* zyzÑÑÑÑ 2(R*− b)/ λ 1 + λ 2(a − R*)/ λ*j o o j z j z j z e 1 1 1 1 1 e − − + − + + − − + j z j j z j z z Å Å Ñ o j z j z j z o b o λ z{ ÅÅÅÅÇk ε* {jk λ z{ jk λ* {ÅÅÅÅjk λ* z{ k ε* {jk λ* z{ÑÑÑÑÖ k o 1− λ o ÅÇÅ o n

( (

) )

( (

( (

) )

) )

( (

The integral in eq 19 needs to be numerically evaluated in order to obtain the electrostatic contribution to the osmotic coefficient. The unknown parameters in the equation are a, b, R* and ϵ*. The parameters a, R*, ϵ*, and b are chosen as temperature-dependent fitting parameters. Note that the value of the parameter, a, the ionic radius, needs to be different from that used in the Debye−Hückel theory. The reason is that since D-H theory ignores the variation in dielectric constant of

) )

| o o o o o o o o o o o o o o } o ÑÉÑo ÑÑo o ÑÑo o ÑÑo o o ÑÑÑo o ÑÑo o ÑÖo ~

(32)

water in the proximity of the ion, it needs to a use a lower value of a than in the model to correct for the stronger electrostatic interaction in the close proximity of the central ion.

4. RESULTS AND DISCUSSIONS 4.1. The Estimation of Interaction Parameters from Regression Analysis. We rewrite the split form of the osmotic coefficient as follows E

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research φ = 1 + φeE + φeNE

(33)

where the first term represents the osmotic coefficient of the ideal solution. The terms φEe , and φNE e , are respectively, the excess electrostatic and the excess nonelectrostatic contributions to the osmotic coefficient. We again rewrite the nonelectrostatic contribution term to osmotic coefficient as follows, φeNE = −

1 [ln(1 − cvh) + c(vh − νvw̅ )] − 1 νmM w

(34) −1

where vh is the molar volume of hydrated salt (dm ·mol ), v̅w is the partial molar volume of free water (dm3·mol−1), and c is the molar concentration of the solution (mol·dm−3). The molar volume of hydrated salt, vh is defined as follows, 3

vh = vs +

hvwb

Figure 2. Plot of excess electrostatic contribution to the osmotic coefficient of aqueous NaCl. Points are from the experiment. Smooth curves are from the present model using parameters given in Table 1. The osmotic coefficient data for NaCl was taken from ref 38.

(35)

vbw

where vs and are, respectively, the molar volume of the bare electrolyte and that of the bound water associated with the electrolyte. The quantity h, represents the primary hydration number of the electrolyte, that is, the number of water molecules tightly bound to the electrolyte. However, vh can be calculated by using slope analysis.1 Since the molar hydrated ion volume, vh is known, eq 34 can be used to estimate φNE e , the nonelectrostatic contribution in the entire concentration range (up to the solubility limit). Now, it is possible to segregate the electrostatic and the nonelectrostatic contributions to the osmotic coefficient in the entire range of electrolyte concentration. Subtracting 1 + φNE e from φ, it is possible to obtain values of φEe , the excess electrostatic contribution (see eq 33). The model developed above in the present paper has four parameters namely a, R*, b, and ϵ*. The nonlinear leastsquares regression method of Levenberg−Marquardt is used to estimate a, R*, b, and ϵ* values. The l2-norm of the residual (∥R∥) is used for judging the quality of the regression. It is defined as ÄÅ ÉÑ1/2 ÅÅ n data i exp model y2 Ñ ÑÑ ÅÅ r − r j z j j j zzz ÑÑÑÑ R = ÅÅÅÅ ∑ jjj zzz ÑÑ σj ÅÅÅ j = 1 jjk { ÑÑÑÖ (36) ÅÇ

The accuracy of the model is illustrated further through the parity plots of the excess osmotic coefficient for aqueous electrolytes in Figure 3. The abscissa represents experimental

Figure 3. Parity plot of the excess osmotic coefficient for aqueous NaCl.

values, and the ordinate represents the values predicted using our model. In Figure 3, the parity is excellent. This model has been tested for other aqueous electrolyte solutions, namely, LiCl, KCl, MgSO4, and MgCl2. Regressed parameters a (nm), R* (nm), b (nm) and ϵ* for these electrolytes are tabulated in Table 2. The excess electrostatic contribution to the osmotic coefficient at various temperatures and the parity plots of these electrolyte solutions are illustrated in Figures 4 to 7. From Tables 1 and 2 and Figure 8a−c, it can be observed that all four parameters namely a, R*, b, and ϵ* have different trends with respect to temperature for different electrolytes. Therefore, as of now, it is out of the scope of the present paper to define a particular function of temperature for all four parameters. However, it is seen from Tables 1 and 2 and Figure 8a−c that, in general, estimated parameters a, R* and ϵ* increase with the increase in temperature. The increase in temperature increases the entropy of water molecules and increases their tendency to leave the shell a, and R*. This tendency is resisted by the electric field. Hence, the temperature dependence is stronger for the electrolyte with greater values of the shell radii ‘a’ and ‘R’ compared to those having lower corresponding values (which translate to lower hydration numbers).

model where rexp (j = 1, ..., n data) respectively represent the j and rj experimental value and the corresponding model prediction of the quantity to be fitted (the experimental data of φEe retrieved from φ, using the method discussed in the aforementioned paragraph) and σj is the standard deviation. Table 1 shows the regression estimates of a(nm), R*(nm), b(nm), and ϵ* for aqueous solutions of NaCl at different

Table 1. Parameters a (nm), R* (nm), b(nm), and ϵ* Used for Aqueous NaCl Solution no.

temp (°C)

a

R*

b

ϵ*

std dev (%)

1 2 3 4

5 25 50 75

0.353 0.367 0.424 0.429

0.455 0.461 0.529 0.535

4.7 4.7 4.7 4.7

25.72 27.04 27.79 28.76

1.11 1.09 1.01 1.72

temperatures in the range from 5 °C to 75 °C. The experimental data for the osmotic coefficient have been obtained from the work of Gibbard et al.38 The model fits fornthe experimental data are shown in Figure 2. The fits are excellent at all temperatures. F

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Table 2. Parameters a (nm), R* (nm), and ϵ* Used for LiCl, KCl, MgSO4, and MgCl2 T (°C)

a

ε*

R*

std dev (%)

a

R*

LiCl 25 40 50 60 70 75 80

0.421

0. 471

27.17

1.97

0.438

0.510

27.78

0.53

0.448

0.550

27.79

0.65

0.412 0.413 0.416 0.418 0.421

0.510 0.515 0.543 0.549 0.569

0.431

0.609

MgSO4 25 40 50 60 70 80

0. 449 0. 459 0.491 0.471 0.472 0. 474

0.500 0.520 0.530 0.555 0.560 0.571

ε*

std dev (%)

21.92 22.91 23.74 26.8 27.83

1.61 1.24 1.98 1.39 1.44

28.81

1.20

KCl

MgCl2 9.99 10.98 11.96 12.92 13.7 14.5

0.86 0.89 0.99 1.14 1.32 1.95

0.499

0.530

6.01

1.74

0.550 0.551 0.562 0.582

0.600 0.605 0.627 0.665

7.29 7.57 8.12 8.96

1.71 1.07 1.82 2.60

Figure 4. (a) Plot of excess electrostatic contribution to the osmotic coefficient for aqueous LiCl; (b) parity plot of the excess osmotic coefficient for aqueous LiCl. The osmotic coefficient data for LiCl was obtained from ref 39.

Figure 5. (a) Plot of excess electrostatic contribution to the osmotic coefficient for aqueous KCl; (b) parity plot of the excess osmotic coefficient for aqueous KCl. The osmotic coefficient data for KCl were obtained from ref 40.

In the former case, the shell a, and R* consists of multiple water layers and the water molecules in the outermost layer are subjected to lower electric field and hence are more prone to leave the hydration shell at a higher temperature. On the

other hand, those ions with very large Pauling radii have perhaps only a single shell of the water molecule, which is strongly bound to the ion and hence the water molecules in that shell are more resistant to detachment at a higher G

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 6. (a) Plot of excess electrostatic contribution to the osmotic coefficient for aqueous MgSO4; (b) parity plot of the excess osmotic coefficient for aqueous MgSO4. The osmotic coefficient data for MgSO4 was obtained from ref 40.

Figure 7. (a) Plot of excess electrostatic contribution to the osmotic coefficient for aqueous MgCl2; (b) parity plot of the excess osmotic coefficient for aqueous MgCl2. The osmotic coefficient data for MgCl2 was obtained from ref 40.

Figure 8. Plot of the models̀ parameters with respect to temperature for aqueous NaCl, LiCl, KCl, MgSO4, and MgCl2.

in the presence of the electric field. In fact, the bulk dielectric constant is the measurement of the orientation of water molecules toward the applied external field: the higher is the value of the dielectric constant, the greater will be the extent of orientation of water molecule toward the applied external field. As temperature increases, the thermal energy of water molecules also increases which in turn increases the vibrational energy and rotational energy of water molecules. Because of the increment in the vibrational energy and rotational energy

temperature. Potassium salts are the example of this category, where the hydration number is insensitive to temperature. Potassium salts have a small hydration number, but the Pauling radius of the potassium ion is 1.33 Å which is larger than most monovalent ions. Surprisingly, the results show that the dielectric constant, ϵ* in the region a < r < R* increases with temperature while the bulk dielectric constant decreases with temperature. Dielectric constant is the measure of the orientation of water molecules H

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research of the water molecule, its degree of freedom also increases which causes the loss in the orientation of water molecules toward the applied external field. Hence the bulk dielectric constant ϵ decreases with increase in temperature. However, in the region a < r < R*, the water molecules are highly orientated toward the central ion because of the presence of a very high electric field of the ion. These water molecules show very less orientation toward the applied external field. Hence they have a low dielectric constant. As temperature increases, the thermal energy of water molecules in the region a < r < R* also increases, which decreases the influence of the electric field generated by the central ion on these water molecules. Hence the influence of the external electric field on them become stronger which results in the higher dielectric constant of these water molecules with increase in temperature. The analysis of the results of the studied electrolytes suggests that our model is limited to the range of concentration of a very dilute solution to a saturated solution and is valid for the temperature range 5 to 75 °C. 4.2. Importance of Hydration Number of Electrolyte Solutions in the Dilute Concentration Region. Among the unknown parameters, a is the hydrated radius which contains the ion (cation/anion) with the primary hydrated shell of water molecules. In the extended Debye−Hückel theory, a is the only interaction parameter. The extended Debye−Hückel theory ignores the nonelectrostatic interaction; hence, the fitted parameter a takes care of both electrostatic and nonelectrostatic interaction. Here, it can be observed that nonelectrostatic interactions are important even in the very low concentration range at which Debye−Hückel theory applies. In Figure 9, we

1 z+z −Aφ I σ(κa) 3

φeDH =

(37)

where I is ionic strength, given by I=

1 c ∑ cizk2 = |z+z−|ν 2 i 2

(38)

κ is the inverse Debye length, given by F I 2εsε0RT

κ=

(40)

and σ(x) is the concentration correction factor, given by σ (x ) =

3 1 (1 + x − − 2 ln(1 + x)) 3 1+x x

(41)

and Aφis the Debye−Hückel constant, given by Aφ =

F 2e 2 8π (εsε0RT )3/2

(42)

Ionic radius a needed in eq 37 is obtained from Stokes and Robinson, 1973.42 In our method, we include both electrostatic and nonelectrostatic contributions to estimate the osmotic coefficient for the dilute solution and hence the total excess osmotic coefficient is given by the following equation φe =

1 1 z+z −Aφ I σ(κa) − [ln(1 − cvh) 3 νmM w + c(vh − νvw̅ )] − 1

(43)

where a is estimated using the procedure described in section4.1 and second term of eq 43 is the nonelectrostatic contribution obtained in the our previous work.1 We have plotted the percentage error between the prediction and the experiments which are estimated using the following equation %error =

(φepredicted − φeexperimental) φeexperimental

100 (44)

The % error associated with the estimation of the excess osmotic coefficient by the two methods are compared in Figure 10 for different electrolytes. We see from this plot that the percent error in our method is less than 5% while that of the Extended Debye−Hückel theory is as high as 30% in some cases. So from this, we conclude that the nonelectrostatic contribution cannot be ignored even at low concentration.

Figure 9. Deviation of excess nonelectrostatic contribution from total excess osmotic coefficient at 25°C

have plotted φNE e /φe versus molality for six different electrolytes in the concentration range 0.001 m to 0.1 m. From Figure 9, it is seen that the nonelectrostatic contribution is a significant part of the total osmotic coefficient. In all cases of the electrolyte solution, the nonelectrostatic contribution becomes more significant with the increase in electrolyte concentration; hence, for a more accurate estimation, both electrostatic and nonelectrostatic should be considered. To prove this point, we have compared the estimated osmotic coefficient of few electrolytes using the extended Debye−Hückel equation and using our method. The details are described below. According to the extended Debye−Hückel theory, the excess electrostatic contribution of the osmotic coefficient is given by the following equation.41

5. CONCLUSION In this work, we have developed the cell model for electrostatic contribution to the osmotic coefficient. From this predictive model, we have regressed the parameters a, R*, ε*, and b. We have made use of a that is, minimum distance approach between the central ion and counterion to estimate the osmotic coefficient of the dilute solution. Here we have also included the nonelectrostatic contribution (which is ignored in D-H theory). After comparison with D-H theory, we have shown that our model is more accurate than the standard D-H theory for NaCl, LiCl, KCl, MgSO4, and MgCl2 solutions. The continuation of this work (correlating the values of a and R* with hydration radius which are estimated in our previous work,1 correlating the I

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 10. (Blue ■) Percentage error predicted using the present method (eq 43), (Red ●) Percentage error estimated using extended Debye− Hückel equation (eq 37). In our method parameters used for electrolytes at 25°C: (a) LiCl [a = 0.421 nm, h = 5.1], (b) LiBr [a = 0.395 nm, h = 5.8], (c) KOH [a = 0.36 nm, h = 5.1 ], (d) HCl [a = 0.46 nm, h = 5.0] .

values of ε* with the Booth model,43 and analysis of temperature dependence of model’s parameters) is under process, and should help us in the future to better understand electrostatic and nonelectrostatic interactions in aqueous electrolyte solutions.



(2) Debye, P.; Hückel, E. The theory of electrolytes. I. Lowering of freezing point and related phenomena. Phys. Z. 1923, 24, 185. (3) March, N. H.; Tosi, M. P. Coulomb Liquids; Academic: London, 1984. (4) Kirkwood, J. G.; Poirier, J. C. The statistical mechanical basis of the Debye-Hückel theory of strong electrolytes. J. Phys. Chem. 1954, 58, 591. (5) Pitzer, K. S. Thermodynamics of electrolytes. I. Theoretical basis and general equations. J. Phys. Chem. 1973, 77, 268. (6) Frank, H. S.; Thompson, P. T. Fluctuations and the limit of validity of the Debye-Hückel theory. J. Chem. Phys. 1959, 31, 1086. (7) Rasaiah, J. C. A view of electrolyte solutions. J. Solution Chem. 1973, 2, 301. (8) Kontogeorgis, G. M.; Folas, G. K. Thermodynamic Models for Industrial Applications. From Classical and Advanced Mixing Rules to Association Theories; J. Wiley and Sons: Chichester: UK, 2010. (9) Lewis, G. N.; Randall, M. In Thermodynamics, 2nd ed.; Pitzer, K.S., Brewer, L. Eds.; McGraw-Hill: New York, 1961. (10) Ciavatta, L. The specific interaction theory in equilibrium analysis. Some empirical rules for estimating interaction coefficients of metal ion complexes. Annal. Chim. (Roma) 1990, 80, 255. (11) Renon, H.; Prausnitz, J. M. Local Compositions in Thermodynamic Excess Functions for Liquid Mixtures. AIChE J. 1968, 14 (1), 135. (12) Bollas, G. M.; Chen, C. C.; Barton, P. I. Refined electrolyteNRTL model: Activity coefficient expressions for application to multielectrolyte systems. AIChE J. 2008, 54, 1608. (13) Chen, C. C.; Mathias, P. M.; Orbey, H. Use of hydration and dissociation chemistries with the electrolyteNRTL model. AIChE J. 1999, 45, 1576.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.8b04986. Detailed derivation of Frank’s model(PDF)



AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. ORCID

Jyoti Sahu: 0000-0003-0485-6298 Notes

The authors declare no competing financial interest.



REFERENCES

(1) (a) Sahu, J.; Juvekar, V. A. Development of a rationale for decoupling osmotic coefficient of electrolytes into electrostatic and nonelectrostatic contributions. Fluid Phase Equilib. 2018, 460, 57. (b) Sahu, J.; Juvekar, V. A. Data on primary hydration characteristics of aqueous electrolytes. Data in Brief 2018, 19, 485. J

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

373 K and from 1 to 6 mol kg−1, and related properties. J. Chem. Eng. Data 1974, 19, 281. (39) Gibbard, H. E.; Scatchard, G. Liquid-Vapor Equilibrium of Aqueous Lithium Chloride, from 25 ̊to 100̊C and from 1.0 to 18.5 Molal, and Related Properties. J. Chem. Eng. Data 1973, 18, 293. (40) Snipes, H. P.; Manly, C.; Ensor, D. D. Heats of dilution of aqueous electrolytes. Temperature dependence. J. Chem. Eng. Data 1975, 20, 287. (41) Newman, J.; Thomas-Alyea, K. E. Electrochemical Systems, 3rd ed.; Prentice-Hall: 1973. (42) Stokes, R. H.; Robinson, R. A. Solvation equilibria in very concentrated electrolyte solutions. J. Solution Chem. 1973, 2, 173. (43) Booth, F. The Dielectric Constant of Water and the Saturation Effect. J. Chem. Phys. 1951, 19, 391.

(14) Chen, C. C.; Song, Y. Generalized electrolyte-NRTL model for mixed-solvent electrolyte systems. AIChE J. 2004, 50, 1928. (15) Chen, C. C.; Britt, H. I.; Boston, J. F.; Evans, L. B. Local composition model for excess Gibbs energy of electrolyte systems. AIChE J. 1982, 28, 588. (16) Chen, C. C.; Evans, L. B. Local composition model for excess Gibbs energy of aqueous electrolyte systems. AIChE J. 1986, 32, 444. (17) Liu, Y.; Harvey, A. H.; Prausnitz, J. M. Thermodynamics of concentrated electrolyte solutions. Chem. Eng. Commun. 1989, 77, 43. (18) Song, W.; Larson, M. A. Activity coefficient model of concentrated electrolyte solutions. AIChE J. 1990, 36, 1896. (19) Kolker, A. R. Thermodynamics modeling of concentrated aqueous electrolyte and nonaqueous systems. Fluid Phase Equilib. 1991, 69, 155. (20) Kolker, A.; Pablo, J. J. D. Thermodynamic modeling of concentrated aqueous electrolyte and nonelectrolyte solutions. AIChE J. 1995, 41, 1563. (21) Nicolaisen, H.; Rasmussen, P.; Sorensen, J. M. Correlation and prediction of mineral solubilities in the reciprocal salt system (Na+, K+) (Cl−, SO42‑)-H2O at 0−100°C. Chem. Eng. Sci. 1993, 48, 3149. (22) Thomsen, K. Modeling electrolyte solutions with the extended universal quasichemical (uniquac) model. Pure Appl. Chem. 2005, 77, 531. (23) Thomsen, K.; Rasmussen, P.; Gani, R. Correlation and prediction of thermal properties and phase behavior for a class of electrolyte systems. Chem. Eng. Sci. 1996, 51, 3675. (24) Lu, X. H.; Maurer, G. Model for describing activity coefficients in mixed electrolyte solutions. AIChE J. 1993, 39, 1527. (25) Fredenslund, A.; Jones, R. L.; Prausnitz, J. M. GroupContribution Estimation of Activity Coefficients in Nonideal Liquid Mixtures. AIChE J. 1975, 21, 1086. (26) Wilson, G. M. Vapor-liquid equilibrium. xi. a new expression for the excess free energy of mixing. J. Am. Chem. Soc. 1964, 86, 127. (27) Hingerl, F. F.; Wagner, T.; Kulik, D. A.; Thomsen, K.; Driesner, T. A new aqueous activity model for geothermal brines in the system Na-K-Ca-Mg-H-Cl-SO4-H2O from 25 to 300°C. Chem. Geol. 2014, 381, 78. (28) Parvaneh, K.; Haghbakhsh, R.; Shariati, A.; Peters, C. J. A crossover UNIQUAC model for critical and noncritical LLE calculations. AIChE J. 2015, 61, 3094. (29) Gebreyohannes, S.; Neely, B. J.; Gasem, K. A. M. Generalized interaction parameter for the modified Nonrandom Two-Liquid (NRTL) activity coefficient model. Ind. Eng. Chem. Res. 2014, 53, 20247. (30) Varela, L. M.; Garcia, M.; Mosquera, V. Exact mean-field theory of ionic solutions: non-Debye screening. Phys. Rep. 2003, 382, 1. (31) Xiao, T.; Song, X. A molecular Debye-Hückel theory and its applications to electrolyte solutions: The size asymmetric case. J. Chem. Phys. 2017, 146, 124118−1. (32) Xiao, T.; Song, X. A molecular Debye-Hückel theory of solvation in polar fluids: An extension of the Born model. J. Chem. Phys. 2017, 147, 214502. (33) Moggia, E. Application of the quasi-random lattice model to rare-earth halide solutions for the computation of their osmotic and mean activity coefficients. J. Rare Earths 2014, 32 (10), 979. (34) Moggia, E. Pseudolattice Theory of electrolyte solutions: Consistency analysis of the Quasi-Random Lattice model at infinite dilution. Fluid Phase Equilib. 2016, 427, 231. (35) Frank, H. S. Local Dielectric Constant and Solute Activity. A Hydration-Association Model for Strong Electrolytes. J. Am. Chem. Soc. 1941, 63, 1789. (36) Bockris, J. O’M., Reddy, A. K. N. Modern Electrochemistry: An Introduction to an interdisciplinary Area, 2nd ed.; Springer-Verlag: New York, 1998. (37) Onsager, L. Theories of Concentrated Electrolytes. Chem. Rev. 1933, 13, 73. (38) Gibbard, H. F., Jr.; Scatchard, G.; Rousseau, R. A.; Creek, J. L. Liquid-vapor equilibrium of aqueous sodium chloride, from 298 to K

DOI: 10.1021/acs.iecr.8b04986 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX