Research Needed: Validity of speciation + Pitzer activity in native engine

[rkingsbury]

Speciation support was recently added to the native modeling engine (in v0.8.0), which up until that point only provided activity coefficients for pre-defined ionic compositions. Hence, speciation and thermodynamics are now combined into one modeling engine.

Although in many cases, users will be primarily interested in one or the other (speciation OR activity coefficients of strong electrolytes), some research is needed to more fully examine the validity of combining these approaches - i.e., what happens when you try to calculate activity coefficients on a solution after speciation?

There are several specific, related items to investigate here:

Pitzer activity coefficients of speciated strong electrolytes

Consider a solution of MgCl2. If I enter this as a simple electrolyte solution, pyEQL's pitzer model implementation will return accurate activity coefficients:

>>> from pyEQL import Solution
>>> s1 = Solution(solutes={'Mg+2':'0.5 mol/kg', 'Cl-': '1 mol/kg'})
>>> s1.get_activity_coefficient('Mg+2')
<Quantity(0.47772752, 'dimensionless')>
>>> s1.ionic_strength
<Quantity(1.5000001, 'mole / kilogram')>

If I speciate the solution, there are many new species besides Mg+2 and Cl-:

s1.equilibrate()
s1.components
{'H2O(aq)': 55.36139419422551, 'Cl[-1]': 0.9166666167185556, 'Mg[+2]': 0.4181428103898319, 'MgCl[+1]': 0.08038100818816704, 'OH[-1]': 1.4955705110947125e-07, 'H[+1]': 1.1847317184369182e-07, 'HCl(aq)': 1.2252437369715798e-08, 'MgOH[+1]': 3.956765514635232e-13, 'O2(aq)': 7.383601288611169e-25, 'HClO(aq)': 1.566134233256325e-27, 'ClO[-1]': 6.446793486822465e-28, 'H2(aq)': 5.759059750780141e-35, 'ClO2[-1]': 0.0, 'ClO3[-1]': 0.0, 'ClO4[-1]': 0.0, 'HClO2(aq)': 0.0}

Now, get_activity_coefficient is programmed to use total salt concentrations rather than individual ion concentrations when invoking the Pitzer model, so the total amount of Mg+2 and Cl- entering the calculation is no different than before speciation. However, the ionic strength is different because many of the new species are monovalent rather than divalent, and hence don't contribute as much. This causes the value of the activity coefficient to shift a little bit.

>>> s1.get_activity_coefficient('Mg+2')
<Quantity(0.473976857, 'dimensionless')>
>>> s1.ionic_strength
<Quantity(1.33835485, 'mole / kilogram')>

Presumably, this result is less accurate, because when Pitzer parameters are fitted to data, I don't think speciation is considered. In other words, the ionic strength used in the expressions always assumes full dissociation of the respective salts. So, how do we deal with this?

Thermodynamic Inconsistency between PHREEQC and pyEQL Pitzer

Speciation calculations in the native modeling engine use PHREEQC's llnl.dat database, because it contains a large number of species (including trace species) and works well at moderate to high salinities (see Lu et al.). However, it does not use the Pitzer model (rather, it uses the B-dot equation; see Table 1).

Meanwhile, the native modeling engine uses the Pitzer model in an attempt to most accurately capture solution thermodynamics at moderate to high concentrations. This obviously results in some inconsistency between the speciation and the activity coefficients.

The question is - what are the practical implications of this inconsistency? Does it introduce meaningful inaccuracies into the properties of interest, or not?