Miyi Li*,Tao Fang

1 School of Chemical Engineering and the Environment,Beijing Institute of Technology,Beijing 100081,China

2 Beijing Institute of Aerospace Testing Technology,Beijing 100074,China

Keywords:Activity coefficient Electrolyte Ion size Poisson–Boltzmann equation

ABSTRACT A rigorous approach is proposed to model the mean ion activity coefficient for strong electrolyte systems using the Poisson–Boltzmann equation.An effective screening radius similar to the Debye decay length is introduced to define the local composition and new boundary conditions for the central ion.The crystallographic ion size is also considered in the activity coefficient expressions derived and non-electrostatic contributions are neglected.The model is presented for aqueous strong electrolytes and compared with the classical Debye–Hückel(DH)limiting law for dilute solutions.The radial distribution function is compared with the DH and Monte Carlo studies.The mean ion activity coefficients are calculated for 1:1 aqueous solutions containing strong electrolytes composed of alkali halides.The individualion activity coefficients and mean ion activity coefficients in mixed solvents are predicted with the new equations.

1.Introduction

The Debye–Hückel(DH)limiting law[1]is a monumental work for electrolyte systems,followed by many equations developed for these systems.Although it gives a good description for physical properties for dilute solutions,it fails at high concentrations for real systems.Another defect is ignoring the size of ions,leading to large deviations from real physical behavior.Fowler and Guggenheim[2]have improved the DH theory with the Güntel berg charging model,with the central ion and ion atmosphere being charged simultaneously.The equation can be applied for electrolyte systems up to 0.1 mol·kg-1,even in mixed solvent systems.Bromley[3]has introduced one more empirical term for ionic strength,which can be applied for solutions at concentrations up to 6 mol·kg-1.Mayer[4]has developed a DH equation with a statistical treatment of ion–ion interactions,based on his theory of clusters of real gases.Stokes and Robinson[5]have proposed the well known ionic hydrate model for aqueous electrolyte systems at high concentrations,assuming that the electrolyte dissociates and reaches ionic solvation equilibria(chemical reaction)in the solution, fitting the ionic hydrate numbers of water molecules and correlating the mean ion activity coefficients of electrolytes from 0.1 to 4 mol·kg-1with experimental data.

In the 1970s,the mean spherical approximation(MSA)[6,7]was derived from the Ornstein–Zernike integral method.It is based on the linear Poisson–Boltzmann equation with an analytical solution for the central ion and surrounded ionic cloud with finite size.The MSA has been applied for the calculation of departures of real systems from ideal ionic solutions,with the solvent as a dielectric continuum,called primitive model.Ion radius and dielectric constant are fitted forosmotic coefficients and activity coefficients of strong electrolyte solutions[8–10].In the correlation,the Pauling diameters are usually chosen for the anion,while cation diameters are set as adjustable parameters for the ionic strength.All these attempts exaggerate the size of smaller cations,which is interpreted as hydration.For larger ions,smaller diameters are obtained,sometimes even smaller than crystallographic diameters,which is difficult to be accepted.Simonin et al.[11,12]have taken into account the concentration-dependent solvent permittivity with the effect on the Helmholtz energy of the system,providing a better description for real electrolyte systems.Besides,the MSA and ion association have been applied together for 1:1 electrolytes in water with adjustable ionic diameters and associate constant[13].The MSA approach is also combined with the modified non-random two-liquid approach[14].It is suitable for the description of mixed-solvent electrolyte systems.

Outhwaite[15,16]has proposed a modified Poisson–Boltzmann(MPB)equation to fit the ion size and exclusive volume of cloud ions.The radial distribution function is represented as a symmetric function when ion species have different valences or ionic radii.Comparisons of the hypernetted chain closure(HNC)integral equation and Monte Carlo(MC)method for systems with different ion sizes[17]indicate that the symmetric Poisson–Boltzmann(PB)expression is viable for the primitive electrolyte model at low concentrations with high quality.Molero et al.[18]have made great efforts to examine the meaning of individual ion activity coefficients using the MPB equation and compared them with the MSA,HNC and MC simulations.It is found that the electric part of activity coefficient for single ion declines monotonically with the increases of ratio of cation to anion ionic size.Kjellander and Mitchell[19]define a long range contribution as a charge distribution like “dressed ion”,called dressed ion theory(DIT),which follows the DH expression.It assumes an effective point charge instead of bare charge of ions,a decay length instead of DH length,and so on.With these renormalized assumptions the entire theory can be cast in a linear PB form.Varela et al.[20]have exploited the dressed ion model to obtain static structure factors of bulk electrolytes or colloid systems and emphasized that the DIT is an exact theory,applicable for all primitive model systems with symmetric and non-symmetric electrolytes.It extends the classical DH theory considerably and the results are valid from low to finite concentrations.This theory is particularly useful in explaining concentrated electrolyte solutions,colloid media,doublelayer,and interfacial phenomena in surfactant physics.On the basis of modern density functional analysis,Abbas and Nordholm[21]have considered long-range electrostatic interactions and short-range ion size effects by a generalized van der Waals analysis,named the corrected Debye–Hückel(CDH).In this theory,all ions are assumed to have the same diameter for symmetric salts,macro-ions or plannersurfaces,with the excluded volume effect vanishing in the linear response.A mean diameter of cation and anion according to the original DH theory is used to fit experimental mean ion activity coefficients and osmotic coefficients up to 1 mol·kg-1.A new approach to consider the ion size of strong electrolyte has been made by Fraenkel[22],using the distance of closest approach between opposite ions in the DH treatment and defining another second distance of closest approach of similar ion,called the Smaller-ion Shell model.It is a DH level modification based on the LPB equation,providing some analytical expressions for real electrolyte systems regarding ion sizes,ion charges,and ion charge asymmetry.Although it is somewhat empirical and complex in two specially defined ionic cases,it can be used directly to correlate experimental activity coefficient data with adjustable parameters instead of numerical commutation such as integral equation methods.

The objective of this work is to develop a new mean ion activity coefficient expression for strong electrolyte systems to describe ion size contribution based on the Poisson–Boltzmann equation.Cations and anions are considered as charged hard spheres of unequal size dissolved in solvents.The ion size is more close to the real crystallographic diameters compared to other theories mentioned above by fitting those as adjustable parameters.The electrostatic potential distribution around the central ion is analyzed by introducing an ion atmosphere assumption similar to that in the DH theory,only the electrostatic forces between ions are considered,and the solvent is viewed as a continuum with dielectric constant.New outer boundary conditions for the PB equation are developed to obtain an analytical solution for electrostatic potential.The effective screening radius is represented as an empirical expression of ionic concentration.The derived mean ion activity coefficient equation is applied for single salt in pure water, fitted with experimental data at 298.15 K.The individual ion activity coefficients and mean ion activity coefficients in mixed solvent electrolyte systems are predicted with our model and parameters.

2.Modeling for the Mean Ion Activity Coefficient

We start from single 1:1 strong electrolytes in pure water,with the salt completely dissociating into ions,which are assumed as hard spheres interacting through Coulomb forces in a dielectric continuum solvent.For simplification,ion size and dielectric constant of solvent do not vary with ionic concentration.Some theoretical relations used here are based on the Debye–Hückel theory[1]and the solutions at low and moderate concentrations are considered.We first recall some of the Debye–Hückel theory to give a specific derivation of our model and systems of alkali halides in water are investigated.

In the bulk thermodynamic relations of the Debye–Hückel theory,the non-ideality of solution is mainly attributed to the ion–ion interactions.The chemical potential for a given ion can be described as a sum of an ideal term and an excess term representing the electrostatic force.The relationship of chemical potential and activity coefficient of ion has been given elsewhere[1,5,24].

Superscripts “ideal”and “electro”refer to the ideal solution and electrostatic energy contribution,respectively,γjis the activity coefficient of ion j in the solution,R(8.314 kJ·kmol-1·K-1)is the gas constant,and T(K)is the general absolute system temperature.

The problem starts from how to determine the contribution of electric energy to the chemical potential.The DH theory calculates the local potential of an ion due to its ionic atmosphere and gives a radial distribution function using the Boltzmann law for ion j in the ionic atmosphere.The same expression for the charge density is used for ion j.

where ρeis the average charge density at distance r from the centralion in molar scale,zjis the valence of ion j,e(1.60218×10-19C)is the fundamental electronic charge,ρj0is the overall average concentration of ion j(mol·cm-3),and ψ indicates the local electrostatic potential outside the central ion.In Eq.(3),the Boltzmann distribution function is used to describe the concentration of ions around the central ion relative to the distance from the central ion.The thermal kinetic energy kT is assumed greater than the electrostatic energy zjeψ(r),similar to that in the DH theory.The term e∑ρj0zjmust be equal to zero,which indicates the electroneutrality of the solution.Using the second term of the Taylor series expansion in Eq.(4),the expression for the charge density becomes

With the charge density ρesimilar to the DH theory,the Poisson–Boltzmann equation can be obtained.

where D is the dielectric constant of solvent given by

εWis the relative dielectric constant of water and εW=78.3,4πε0=1.11265 × 10-10C2·N-1·m-2,and the Debye lengthis given by

Substituting Eqs.(5)and(8)into Eq.(6),an analytical solution of the PB equation is obtained.

where parameters A and B can be determined from the boundary conditions.In the DH theory,electrostatic potential ψ approaches zero at infinity distance,so that parameter B is zero.However,this boundary condition is not appropriate for moderate and high concentrations,since the zero potential boundary would move from the in finite distance to the central ion as concentration increases.It seems that the contribution of the last term in Eq.(9)should be introduced.A specific analysis will be performed to obtain appropriate boundary conditions for our modeling.

In electrostatic problems,it is important to determine the electric potential and electric field for a certain charge distribution.We start from the simplest case with a cation+q at(-1,0)and an anion-q at(+1,0)in the same size as shown in Fig.1(a).The electrostatic field can be described by using the equal electric potential surface.When the middle point O between+q and-q is used as the origin,the potential at an arbitrary point P in the field is in accordance with the superposition principle of electric potential[23].

Fig.1.The scheme of electrostatic potential.(a)the simplest case;(b)the second case.

Another case is given by a system consisting of a few charges with anions and cations in a local composition,as shown in Fig.1(b).With the assumptions above,ions with the same charge may not have the density of Eq.(3).By viewing the local composition as a sphere with origin O,the charge function is presented using a charge shell of q with differential elements d r.

At Point P,the electric potential can be obtained using Eq.(10)by accumulating all the contribution of ions in this sphere.The whole system is separated into many realms with this composition assumption,as shown in Fig.2.We define an effective screening radius τ to present this local realm and make it a different meaning from the DH theory.Two boundary conditions are used to determine τ.

Fig.2.The scheme of effective screening radius related to ionic size.

It indicates the superposition principle of electric potential and Eq.(12)presents a zero potential boundary.The total charge outside the central ion must equal its opposite charge in the local area and Eq.(13)indicates the local electric neutrality.It is easily deduced from the Gauss law that the net flux of electric field through the boundary surface at r=τ is zero,which indicates a screening mechanism.Thus parameters A and B become

In these equations,the Debye lengthis well defined and can be obtained in literature[1,5,24,25].Unlike other models,the Pauling radii σjof anions and cations are directly used as strong physical properties here.The new defined parameter τ will be discussed in Section 3.To proceed further,one must separate the contribution to ψ(r)due to the ionic atmosphere from the contribution by the central ion in the absence of other ions,which is the so-called self-atmosphere potential:

When the ionic atmosphere potential contribution φ is estimated at r=σjon the surface of central ion j,we obtain

The remaining work involves relating the energy work in forming the ionic atmosphere to activity coefficient γj.When the departure from ideality is attributed to the effects of ion–ion interactions in Eq.(1),the energy work is estimated as a difference in which the electrostatic potential energy(Wje)of the central ion in the presence of all the other ions with the in finite dilute solution as the reference state.

where k=1.38 × 10-23J·K-1,NA=6.022 × 1023,I(mol·kg-1)indicates the ionic strength,κDHis the Debye decay length in molal unit,and dw(kg·m-3)is the pure solvent density at system temperature and pressure.The mean ion activity coefficient is defined as

where ν+and ν-denote the stoichiometric coefficients of cation and anion,respectively,ν+= ν-=1;z+and z-denote the valence of cation and anion,respectively,z+=+1 and z-=-1 for alkali halides.By substituting the Pauling radius and potential Eq.(17)into Eq.(22),we obtain

Using Eq.(23)with Eqs.(9),(14),(15),and(17),the mean ion activity coefficients can be calculated.

3.Results and Discussion

The Debye–Hückel model provides a rigorous theory for diluted electrolyte solutions.The new mean ion activity coefficient expression for dilute,moderate and higher concentrations obtained in this work is checked with the DH theory in dilute solutions.

3.1.Comparison with the DH and Monte Carlo theory

In the low concentration range with I＜ 0.1 mol·kg-1,the limitation can be obtained withτ+?σ+,τ-?σ-,andκDHτ-≈κDHτ+?1,where σ+andσ-are the crystallographic size of positive and negative ions,respectively.With these relations,Eqs.(14)and(15)become

Using Eq.(9),the electrostatic potential is simplified to

where the electric potential expression outside the central ion is similar to that of the DH theory,and parameterσjindicates the Pauling radius of centralion instead of the nearest contact distance of the opposite ions in the DH theory,which is usually set to the value of 0.400 or 0.425 nm[2,38].Although it is greater than the Pauling radius of ions,its influence on the activity coefficients of ions is neglected for dilute solutions in the DH limiting law[1].In this work,cations Li+,Na+,K+,and Rb+and anions Cl-,Br-,and I-are chosen to verify the expressions presented here.Pure water is used as solvent.1:1 electrolytes composed of these ions in aqueous binary systems are used at 298.15 K and atmospheric pressure.The Pauling radii are taken from Ref.[25]and listed in Table 1.For alkali metals and halogen elements,the radius of ion increases with the index of element.Alkali metal ion radii are smaller while halogen radii are larger than atom radii.

Parameter τjis the effective screening radius proposed in this work.It is a hypothetical variable to define the zero potential of screening local components.No experimental data are available for comparison,so it is expressed as an empirical expression related to the DH decay length and the Pauling radius:

Table 1 The crystallographic size of ions[25]

where σjis the central ion radius and σiindicates the unlike charge ion radius for 1:1 electrolyte systems.With the expression of mean ion activity coefficient,effective screening radius,and parameters of the Pauling radii,12 aqueous systems are investigated.Three adjustable parameters a,b,and c in Eq.(27)are listed in Table 2.

Another approach to check our model is the radial distribution function

where gji(r)is the radial distribution function with ion j as the center and i indicates unlike and like charge ions in the effectivescreening radius.ψ(r)is calculated with Eqs.(9),(14),(15),and(27).The experimental results of the primitive model for an electrolyte solution are provided by the MC study,carried out by Card and Valleau[28].

Table 2 Parameters a,b,and c(d=1/3)

Fig.3.Comparison of the radial distribution functions.

Fig.3 compares our results with those from the Debye–Hückel theory and MC simulations.At lower concentration c=0.1038 mol·L-1,the DH theory,MC simulations and the expression in this work agree reasonably well,indicating that our model obeys the DH theory for dilute and low concentration solutions.At higher concentration c=1.968 mol·L-1,g--and g++match the DH and MC values better,but g+-and g-+show distinctly smaller values than those with the DH theory,while MC data present larger values.It is because the MC and DH theories are carried out from the nearest contact distance 0.425 nm while this work uses gji(r)that is from the radius of central ion j.The boundary conditions are different with the DH theory.In dilute solutions,the effective screening radius has a large value and is in agreement with the DH law,while at high concentrations,the tendency is in agreement.Using the parameters correlated from experimental results in Table 2,the effective screening radii for common cation and anion are investigated using Eq.(27).The DH decay length is listed against the molal concentration of solutions in Table 3.

For halogen anions,the effective screening radius decreases as the central ion size increases at dilute to high concentrations.For alike metal cations,the effective screening radii increase as the central ion size decreases in low concentration solutions,which can be explained with the ionic hydration theory by Robinson and Stokes[5]that smaller cations always have a larger hydration layer.At concentrations higher than 0.5 mol·kg-1,the ionic hydration layer may be weakened and the cation is smaller with smaller effective screening radius.For all concentrations considered,the DH length is always smaller than the effective screening length.Here,the effective screening radius indicates the distance of zero electric potential from the central ion.It has a considerable value compared with the in finite distance of the DH theory with reference zero potential.

To simplify the correlation,parameter d for all the ions is set to 1/3,parameters b and c for larger ions are set to constant 1.50 and 2.00 as shown in Table 2.Parameter a is correlated with the experimental mean ion activity coefficient data in literature[29].For cations Li+and Na+with smaller ionic size,parameters a,b,and c are correlated separately.The calculated mean ion activity coefficients of four salts with common anion Cl-are shown in Fig.4 together with the experimental values.The mean ion activity coefficients increase in the order RbCl＜KCl＜NaCl＜LiCl.The calculated results are in good agreement with the experimental data in the dilute and saturated concentration ranges.It is found that smaller cation size and effective screening radius lead to greater mean ion activity coefficients at the same molality concentration.

3.2.Application to 1:1 aqueous electrolyte systems

The results for some alkali salts with a common anion are shown in Fig.5,and the average absolute relative deviations(AARDs)are listed in Table 4.

where Ndataindicates the number of experimental data,and superscripts“calc”and “exp”refer to calculated and experimental values,respectively.The mean ion activity coefficients increase in the order MCl＜MBr＜MI,where“M”indicates alkali metals Li,Na,K,and Rb.Calculated results agree with the measurements with high accuracy at dilute and moderate concentrations,as well as with the saturated NaCl+H2O and KCl+H2O systems at 298.15 K.Thus this model gives good prediction.With the increase of concentration and ionic size,calculated results are clearly off the experimental results,because larger ions and higher ionic strength have a greater influence on molar volumes of solutions[30],which is neglected in this work.Moreover,only the electrostatic force is considered in the modeling,while solvent–ion and solvent–solvent interactions are ignored.The dielectric parameter of the solvent is influenced by high concentrations[31],which is set as constant here.These simplifications may be responsible for the discrepancy from real chemical systems at relatively high concentrations.

Table 3 The effective screening radius of central ions with common cation and anion

Fig.4.Mean ion activity coefficients of 1:1 electrolytes in water with common anion.Line:our model;symbol:data in literature[29].

3.3.Correlation of the individual ion activity coefficient

Vera et al.[32,33]have demonstrated that the activity of individual ions is a measurable property.Reliable values can be obtained in electrochemical cells with a half-cell ion selective electrode and a standard single junction reference electrode immersed in the sample solution.Some experimental data of individual ion activity coefficients at 298.15 K have been published.In this model,the mean ion activity is composed of the individualion activity coefficient in Eq.(23),expressed with Eq.(18).Fig.6 shows that the individual ion activity coefficients from our model are in agreement with measurements.For system NaCl+H2O,the mean ion activity coefficient values of Vera et al.are systematically 0.5%higher than those of Hamer and Wu[29].Our results agree very well with Vera et al.'s values for NaCl and NaBr aqueous systems with AARD less than 2%,with a greater deviation for individual ion Cl-in system KCl+H2O.It can be figured out that the individual ion activity coefficients are related to the size of other ions.The values of Na+in system NaBr+H2O are obviously higher than that in system NaCl+H2O,and the values of Cl-in system KCl+H2O are also higher than that in system NaCl+H2O.Our model can present this tendency of in dividualion activity coefficients.Because of limited experimental data for individual ion activity coefficients,a further comparison is not available.

Fig.5.Mean ion activity coefficients of 1:1 electrolytes in water with common cation.Line:this work;symbol:data in literature[29].

Table 4 Average absolute relative deviation of this model

Fig.6.The individual ion activity coefficients of 1:1 electrolytes in water.Line:our model;symbol:data in literature[28,32].

3.4.Prediction for single and mixed solvent electrolyte systems

The mean ion activity coefficients in water–alcohol mixtures are also considered.Available data of NaCl,NaBr and KCl[34–37]in water–methanol and water–ethanol mixtures are used to check the predictability of this model.By assuming the solvent as a continuum medium for the electrolyte systems,the mixing rules for mixed solvent are

The correlated results are compared with the experimental data in Fig.7.The prediction is excellent for the mixed solvent system at low concentrations with larger water mass fraction.The mean ion activity coefficients decrease greatly as the salt free water composition decreases.Larger deviations are observed for the system with higher organic mass fraction.For the arbitrary parameter in Eq.(25),larger deviations may be caused by:(1)the mixing rule for mixed solvent and ionic concentration independent of dielectric constant;(2)parameters a,b,and c for aqueous systems;(3)small values of mean ion activity coefficients;and(4)measurement uncertainties.Predicted results match the data well and show appropriate tendency.Overall,the AARD of our results is less than 4.6%.

4.Conclusions

A new activity coefficient model for aqueous strong electrolyte systems is presented where ion size parameters and new boundary conditions are used to obtain an analytical solution of the PB equation.A new electrostatic potential equation is used and the radial distribution functions match the MC and DH theories at low concentrations.

With the assumption for effective screening radius,activity coefficient expressions are given.Calculated results agree with the experimental data for dilute,moderate,and even high concentrations(＞6 mol·kg-1)with acceptable deviations.The mean ion activity coefficients of alkali halides with common anion decrease with increasing size of cation,and the values with common cation increase with the size of anion.

The individual ion activity coefficients are represented with our model for aqueous electrolyte systems at 298.15 K.The activity coefficients of common cation increase with anion size and those of common anion increase with cation size.This tendency can be well reproduced with the derived expressions.Water–alcohol electrolyte systems are predicted and a good tendency is obtained compared with the experimental data.It proves that the model can be used for prediction of electrolyte systems.

Acknowledgments

The authors would like to thank the DDBST GmbH(Oldenburg,Germany)for providing the Dortmund Data Bank for the modeling.The corresponding author(Miyi Li)is grateful to Prof.Gmehling for the direction of this topic and providing the financial support for his stay in the University of Oldenburg,Germany.

Fig.7.Prediction of mean ion activity coefficients in mixed solvent electrolyte systems.(a)NaCl+H2O+MeOH[33];(b)NaCl+H2O+EtOH[34];(c)NaBr+H2O+EtOH[35];(d)KCl+H2O+MeOH[36].Line:the model in this work;symbol:data in literature.