Prediction of activity coefficients in water-methanol mixtures using a generalized Debye-Huckel model

We propose a generalized Debye-Huckel model from Poisson-Fermi theory to predict the mean activity coefficient of electrolytes in water-methanol mixtures with arbitrary percentage of methanol from 0 to 100%. The model applies to any number of ionic species and accounts for both short and long ion-ion, ion-water, ion-methanol, and water-methanol interactions, the size effect of all particles, and the dielectric effect of mixed-solvent solutions. We also present a numerical algorithm with mathematical and physical details for using the model to fit or predict experimental data. The model has only 3 empirical parameters to fit the experimental data of NaF, NaCl, and NaBr, for example, in pure-water solutions. It then uses another 3 parameters to predict the activities of these salts in mixed-solvent solutions for any percentage of methanol. Values of these parameters show mathematical or physical meaning of ionic activities under variable mixing condition and salt concentration. The algorithm can automatically determine optimal values for the 3 fitting parameters without any manual adjustments.


Introduction
Water and alcohol are ubiquitous and complicated liquids [1,2,3]. With salts, they pose many challenges to thermodynamic modeling for a great variety of applications in a wide range of systems and conditions [4,5,6,7,8,9]. One of the major difficulties for numerous models [9] is to deal with the combinatorial explosion of empirical parameters up to tens of thousands [4] to calculate activity coefficients of electrolyte solutions with different compositions at variable temperature and pressure. Even worse, many parameters do not have physical meaning [10] or offer mathematical hint [4,11] to use.
Generalized Debye-Hückel (DH) models [11,12,13,14] from a Poisson-Fermi (PF) theory [15,16,17,18] developed recently can ease some of these difficulties. These models use only 3 empirical parameters having both physical and mathematical properties to well fit experimental activity data of multicomponent electrolyte solutions in a range of concentrations, temperatures, and pressures. It is shown in [11] that the generalized DH model differs much from Hückel's model [19] (and numerous DH models extended from it since 1925) as their approximations of Born solvation energies are inverse of each other in terms of parameters, which explains why extended DH models need more parameters generally without physical meaning. The PF theory treats ions and water (solvent) molecules of any volume and shape with interstitial voids, and accounts for polarization of water, both short and long ranges of ion-ion and ion-water interactions and correlations, and the non-uniform dielectric response (permittivity) of electrolyte solutions.
We propose here a generalized DH model to predict mean activity coefficients of electrolytes in water-methanol mixtures with any percentage (mole fraction) of methanol x in [0, 1]. The model first uses 3 parameters α H 2 O j for j = 1, 2, 3 to best fit the experimental activity data of NaF [20], NaCl [21], and NaBr [22], for example, in pure-water solutions (i.e., x = 0) using the method of least squares. It then uses another 3 parameters ∆α j to predict the activities in mixed solutions for any arbitrary x 0. The parameters α H 2 O j define a factor function θ(I) of the variable ionic strength I of the solution that in turn modifies the experimental Born radius R 0 i [23,24,25] of an ion i in pure solvent (i.e., I = 0) to an unknown Born radius R Born i = θ(I)R 0 i for any I 0. The other ∆α j are defined by α j + x∆α j for any mixing x. Therefore, all these 6 parameters have clear physical meaning in terms of Born energy. We also provide numerical evidence that their values offer novel hints for future studies on different solutions for which a numerical algorithm is given to show how to implement the model with details. This predictive model and the algorithm can be straightforwardly applied to electrolyte solutions with multi-valent ions, mixed salts, variable temperature, and variable pressure [11].

Theory and algorithm
For K species of ions in water-methanol mixed solvents, the entropy model proposed in [17,18] treats ions, water (denoted by K + 1), and methanol (K + 2) as nonuniform spheres with interstitial voids (K + 3). The total volume V of the system is where v i is the volume of each i th species particle, N i is the total number of all i th species particles, and V K+3 denotes the total volume of all the voids. In bulk solutions, we have the bulk concentrations C B i = N i V and the bulk volume fraction of voids i . If the system is spatially inhomogeneous with variable electric or steric fields, as in realistic systems, the constant C B i then changes to a function C i (r) so does Γ B to a void volume func- in the solvent domain Ω s . It is shown in [18] that the distribution (concentration) of particles is of Fermi-like type, i.e., C i (r) < 1 v i for any arbitrary (or even infinite) electric potential φ(r) at any r ∈ Ω s for all i = 1, · · · , K + 2, where β i = q i /k B T with q i being the charge on species i particles and q K+1 = q K+2 = 0, k B is the Boltzmann constant, T is an absolute temperature, and v 0 is a unit volume. The steric potential S (r) is an entropic measure of crowding or emptiness of particles at r [15,17,18]. If φ(r) = 0, then Γ(r) = Γ B and hence S (r) = 0. The factor v i /v 0 in (2) shows that the steric energy −v i v 0 S (r)k B T of a type i particle at r depends not only on the steric potential S (r) but also on its volume v i similar to the electric energy β i φ(r)k B T depending on both φ(r) and q i [18].
The activity coefficient γ i of an ion of species i in electrolyte solutions describes the deviation of chemical potential of the ion from ideality (γ i = 1). The excess chemical potential µ ex i = k B T ln γ i can be calculated by [12] where φ(r) is the electric potential generated by the ion i in the system domain Ω = Ω i ∪ Ω sh ∪ Ω s shown in Fig. 1, Ω i is the spherical domain occupied by the ion, Ω sh is the solvation shell domain of the ion, Ω s is the rest of solvent domain, 0 denotes the center of the ion, and φ 0 (r) is a potential function when the solvent domain is ideal, i.e., C B j = 0 for all j. The potential function φ(r) can be found by solving the PF equation [11,12,13,14] and the Laplace equation where s = 0 x is the permittivity of a bulk mixed solvent in Ω sh ∪ Ω s , 0 is the vacuum permittivity, x is the dielectric constant of the bulk solvent, l c = √ l B l D /48 is a density-density correlation length [26], and l B and l D are the Bjerrum and Debye lengths, respectively. The dielectric operator s (l 2 c ∇ 2 − 1) yields the permittivity of the electrolyte solution and the polarization of the mixed solvent as functions of r [18].
The same derivation steps developed in [11] apply to Eqs. (4) and (5) for an approximate and analytical solution, i.e., (i) linearization of the fourth-order PF (same as Poisson-Bikerman used in [11]) Eq. (4) for a general binary (K = 2) electrolyte C z 2 A z 1 in the mixed solvent with the valences of the cation C z 1 + and anion A z 2 − being z 1 and z 2 , respectively, (ii) determination of global solutions of the linear PF and Laplace equations in the spherical domain in Fig. 1, and (iii) determination of a unique solution of these two equations with the same set of the interface and boundary conditions proposed in [11]. The analytical solution of Eqs. (4) and (5) is in Ω s , where r = |r|, andC B 4 denote respectively the maximal bulk concentrations of water and methanol considered in this work, and the Debye l x D  and correlation l x c lengths have been generalized to include all particle volumes as shown in Λ x . All formulas are in the same form as those in [11] generalized to include mixtures.
Since the solvation free energy of an ion i varies with salt concentrations, the Born energy in a pure mixed solvent (C B j = 0) should be modified to vary with C B j ≥ 0 for j = 1, 2 [11]. Here, the constant Born radii and in pure water and pure methanol are obtained from experimental hydration Helmholtz free energies ∆H 0 i,H 2 O and ∆H 0 i,MeOH [23,24,25], respectively, as given in Table 1 with other physical values. The effective Born radius R Born i,x in (6) is thus a function of C B j that can be modeled by the simple formula [12] R Born i, where I = 1 2 Σ j C B j z 2 j is the ionic strength of the solution, I = I/M is a dimensionless ionic strength, M is molarity, and α x 1 , α x 2 , and α x 3 are parameters for modifying the Born radius R 0 i,x to fit experimental activity coefficients that change with I.
From (3) and (6), we thus obtain a generalized activity coefficient for each ion i in the mixed electrolyte solution and the mean activity coefficient where + and − denote C z 1 + and A z 2 − , respectively. We use three parameters α H 2 O j for j = 1, 2, 3 to fit the experimental activities of NaF [20], NaCl [21], and NaBr [22], for example, in pure-water solutions. We then use another three to predict the activities of these salts in mixed-solvent solutions for any x in [0, 1]. This equation is derived from α Since the steric potential takes particle volumes and voids into account, the shell volume V sh of the shell domain Ω sh can be determined by the steric potential where O ± x is the occupation (coordination [29,30]) number of solvent molecules in Ω sh , v x = (1 − x)v 3 + xv 4 , and C B x = (1 − x)C B 3 + xC B 4 . The shell radius R sh i,x is thus determined by O ± x . The following algorithm summarizes the proposed model with more details in numerical methods and implementation, where Steps 1 -5 are for fitting and 6 -8 for prediction. Figure 2 shows the mean activity coefficients of NaF, NaCl, NaBr fitted (red curves) and predicted (blue curves) by (16) with α H 2 O j and ∆α j in (17) to experimental data (symbols) from [20,21,22] at x = 0 (in pure water) and x = 0.2, 0.4, 0.6, 0.8, 1 (in mixture or pure methanol), respectively. The values of α H 2 O j and ∆α j for j = 1, 2, 3 are given in Table 2  and |∆α 1 | > |∆α 2 | > |∆α 3 |, which implies the order α x 1 > α x 2 > α x 3 as well from (17). This numerical order gives mathematical hints to these parameters for further use of our model in different conditions or for other electrolyte systems. Figure 3 shows that the factor θ in the effective Born radius θR 0
percentage x (= 0 for red solid curves and 0.2, 0.4, 0.6, 0.8, 1 for others) of methanol in mixtures. Therefore, the parameters α x j have physical meaning in Born energy and their values from Table 2 are in accord with experimental solvation energies of these ions in water-methanol mixtures [32,33], i.e., the Born energy of these anions in H 2 O-MeOH mixture is larger than in pure H 2 O in the same conditions, see e.g. Fig. 3 in [33].
We make some more remarks on the results, model, and algorithm as follows: (i) Our model can fit any set of experimental activity data points as shown by the red curves in Fig. 2 with only 3 parameters (α H 2 O j , j = 1, 2, 3) for which their values are determined automatically (not manually) by the algorithm. It nevertheless requires many precise physical values from experimental sources as shown in Table 1 and the algorithm.
(ii) We manually adjusted the values of ∆α j in Table 2 for prediction. In addition to the properties of significant order and Born energy, these values are not arbitrary but can be verified with experimental data. For example, α  Table 2, and 0.0501 is from fitting MeOH data (not shown). This implies that the value 0.027 is reasonable and verifiable.
(iii) As shown and discussed in [11,12,13,14], the principal determinant of ionic activities by (15) is the effective Born radius R Born i,x (I) (or the ion domain Ω i in Fig. 1) due to the singular charge q i δ(r − 0) of the ion, which is infinite at 0 and thus critically affects γ x ± . The secondary part in (15) is R sh i,x that defines the shell volume of Ω sh (18) in which the electric potential function φ(r) in (6) decreases exponentially to the solvent domain Ω s , see e.g. Fig. 6 in [11]. The tertiary part is φ(r) in Ω s that is derived by the PF Eq. (4). Therefore, our model describes the effects of ion and solvent sizes, ion-ion and ionsolvent interactions, and solution permittivity in this order of 3 determinants and of 3 subdomains.
(iv) We chose and fixed the coordination number O ± x = 18 in (18) for simplicity to produce all the above results because of Figure 2: Mean activity coefficients of NaF, NaCl, NaBr fitted (red curves) and predicted (blue curves) by (16) with α H 2 O j and ∆α j in (17) to experimental data (symbols) from [20,21,22] at x = 0 (in pure water) and x = 0.2, 0.4, 0.6, 0.8, 1 (in mixture or pure methanol), respectively. its secondary effect on γ x ± . Its value can be chosen more precisely from experimental sources for different ions in different conditions [29,30]. It can also be changed to a more specific  (18) and Newton() in the algorithm more complicated for implementation.
(v) The activity Eq. (16) is derived from the first principle volume Eq. (1) which is a foundational proposition of our theory that defines the steric potential in (2) in terms of voids. The steric potential is thus a mean-field summary of all kinds of interactions between any pair of particles in a system such as Coulomb (long range), van der Waals (short), or Lennard-Jones (short) interactions [18] that produce the voids and hence the pressure of the system. Therefore, Eq. (16) does not need any mixing and combining rules (yielding more empirical parameters [9]) for these short-range interactions, and can apply to systems under variable temperature or pressure condition [11]. Furthermore, Eq. (16) accounts for variable permittivity of electrolyte solutions with a dielectric function (r) = 1 in Ω i [11], = x in Ω sh , and = H 2 O C 3 (r) in Ω s [11].
(vi) The model is usually expressed in dimensionless form in implementation with the scaling factors s 1 = e/(k B T ) and s 2 = Å 2 e 2 /(k B T ) for the potential φ(r) and concentration C i (r) variables, i.e., s 1 φ and s 2 C i are dimensionless. Our code of the algorithm is accessible at https://github.com/JinnAIGroup for verification and further development.

Conclusion
We proposed a generalized Debye-Hückel model for calculating and studying the activity of electrolytes in watermethanol mixtures for any number of salt types with arbitrary percentage (mole fraction) of methanol. The model is based on the Poisson-Fermi theory that accounts for the effects of (i) non-uniform sizes of ions and solvents, (ii) short and long interactions between ion and solvent or different ions or different solvents by mean-field steric and electric potentials, and (iii) non-uniform and size-dependent permittivity of the mixed solution.
We also proposed an algorithm to implement the model that can automatically and well fit any set of experimental activity coefficients with corresponding salt concentrations using only 3 empirical parameters that show clear physical meaning in terms of Born energy and the significant order of their values for verification and numerical hints for further applications to other electrolyte systems. Based on these parameters, the algorithm can also predict the activity of mixtures using another 3 parameters for any mole fraction of a solvent to another solvent. Again, the later 3 parameters have the same physical meaning and significant order, and are verifiable with experimental data. Our model and algorithm with the same parameters can straightforwardly apply to other electrolyte systems for both fitting and prediction under different conditions such as temperature or pressure.