Ionic strength reported and manually calculated are not consistent

[klk]

Dear David and Collaborators,

When using PHREEQC to perform speciation calculations for aqueous solutions, I noticed that the reported ionic strength does not exactly match the manually calculated value using the molalities of the species and their charges. The following example demonstrates this more significantly (a 4.6% mismatch, with reported value being 4.198 and expected value 4.400).

While this mismatch may be considered low for practical reasons, it significantly affects the computed activities and activity coefficients of some species (e.g. la_Na+ should be 0.0258389, but PHREEQC reports 0.0190442 and therefore a mismatch of 35.7%).

PHREEQC version: 3.7.3
Database: lnll.dat

Code: [Select]

SOLUTION
    temperature 60.0
    pressure 100.0
    units mol/kgw
    pH 7.0 charge
    B  0.10
    Ba 1.20
    Br 0.20
    C  4.00
    Ca 0.50
    Cl 2.00
    Fe 0.70
    K  0.80
    Li 0.60
    Mg 0.40
    Mn 0.60
    Na 2.00
    S  0.50
    Si 0.10
    Sr 1.40
SELECTED_OUTPUT
    -file complex-solution.out
    -high_precision true
    -ionic_strength
    -molalities OH- H+ H2O BH4- BO2- NaB(OH)4 BaB(OH)4+ MgB(OH)4+ B(OH)3 CaB(OH)4+ B2O(OH)5- BaCO3 Ba+2 BaCl+ BaB(OH)4+ BaOH+ Br3- Br- NaBr KBr Br2 BrO- HBrO BrO3- BrO4- C2H4 C2H6 CH4 CO SrCO3 BaCO3 MnCO3 CaCO3 MgCO3 CO3-2 HCO3- NaHCO3 NaCO3- MgHCO3+ CaHCO3+ MnHCO3+ CO2 FeCO3 FeHCO3+ FeCO3+ CaCO3 Ca+2 CaSO4 CaHCO3+ CaB(OH)4+ CaCl+ CaCl2 CaOH+ Cl- NaCl KCl SrCl+ BaCl+ LiCl MgCl+ MnCl+ CaCl+ CaCl2 MnCl3- FeCl+ FeCl4-2 HCl FeCl2 FeCl2+ FeCl+2 FeCl4- ClO- HClO ClO2- HClO2 ClO3- ClO4- FeCO3 FeHCO3+ Fe+2 FeSO4 FeOH+ FeCl+ FeCl4-2 Fe(OH)2 FeCl2 Fe(OH)3- Fe(OH)4-2 Fe(OH)4- Fe(OH)3 Fe(OH)2+ FeOH+2 FeCO3+ Fe3(OH)4+5 Fe2(OH)2+4 FeCl2+ Fe+3 FeCl+2 FeSO4+ Fe(SO4)2- FeCl4- H2 K+ KSO4- KCl KBr KOH KHSO4 Li+ LiSO4- LiCl LiOH MgCO3 MgSO4 Mg+2 MgHCO3+ MgB(OH)4+ MgCl+ Mg4(OH)4+4 MnCO3 Mn2(OH)3+ MnCl+ MnSO4 Mn+2 MnHCO3+ MnCl3- MnOH+ Mn2OH+3 Mn(OH)2 Mn(OH)3- Mn(OH)4-2 Mn+3 MnO4-2 MnO4- Na+ NaCl NaSO4- NaHSiO3 NaHCO3 NaCO3- NaB(OH)4 NaBr NaOH O2 HS- S-2 H2S S2-2 S3-2 S4-2 S5-2 S2O3-2 HS2O3- S2O4-2 SO3-2 HSO3- SO2 H2SO3 S2O6-2 S3O6-2 S4O6-2 S5O6-2 S2O5-2 SO4-2 NaSO4- KSO4- SrSO4 MgSO4 LiSO4- CaSO4 MnSO4 FeSO4 HSO4- KHSO4 FeSO4+ Fe(SO4)2- H2SO4 S2O8-2 HSO5- NaHSiO3 HSiO3- SiO2 H4(H2SiO4)4-4 H6(H2SiO4)4-2 H2SiO4-2 SrCO3 Sr+2 SrSO4 SrCl+ SrOH+
    -activities OH- H+ H2O BH4- BO2- NaB(OH)4 BaB(OH)4+ MgB(OH)4+ B(OH)3 CaB(OH)4+ B2O(OH)5- BaCO3 Ba+2 BaCl+ BaB(OH)4+ BaOH+ Br3- Br- NaBr KBr Br2 BrO- HBrO BrO3- BrO4- C2H4 C2H6 CH4 CO SrCO3 BaCO3 MnCO3 CaCO3 MgCO3 CO3-2 HCO3- NaHCO3 NaCO3- MgHCO3+ CaHCO3+ MnHCO3+ CO2 FeCO3 FeHCO3+ FeCO3+ CaCO3 Ca+2 CaSO4 CaHCO3+ CaB(OH)4+ CaCl+ CaCl2 CaOH+ Cl- NaCl KCl SrCl+ BaCl+ LiCl MgCl+ MnCl+ CaCl+ CaCl2 MnCl3- FeCl+ FeCl4-2 HCl FeCl2 FeCl2+ FeCl+2 FeCl4- ClO- HClO ClO2- HClO2 ClO3- ClO4- FeCO3 FeHCO3+ Fe+2 FeSO4 FeOH+ FeCl+ FeCl4-2 Fe(OH)2 FeCl2 Fe(OH)3- Fe(OH)4-2 Fe(OH)4- Fe(OH)3 Fe(OH)2+ FeOH+2 FeCO3+ Fe3(OH)4+5 Fe2(OH)2+4 FeCl2+ Fe+3 FeCl+2 FeSO4+ Fe(SO4)2- FeCl4- H2 K+ KSO4- KCl KBr KOH KHSO4 Li+ LiSO4- LiCl LiOH MgCO3 MgSO4 Mg+2 MgHCO3+ MgB(OH)4+ MgCl+ Mg4(OH)4+4 MnCO3 Mn2(OH)3+ MnCl+ MnSO4 Mn+2 MnHCO3+ MnCl3- MnOH+ Mn2OH+3 Mn(OH)2 Mn(OH)3- Mn(OH)4-2 Mn+3 MnO4-2 MnO4- Na+ NaCl NaSO4- NaHSiO3 NaHCO3 NaCO3- NaB(OH)4 NaBr NaOH O2 HS- S-2 H2S S2-2 S3-2 S4-2 S5-2 S2O3-2 HS2O3- S2O4-2 SO3-2 HSO3- SO2 H2SO3 S2O6-2 S3O6-2 S4O6-2 S5O6-2 S2O5-2 SO4-2 NaSO4- KSO4- SrSO4 MgSO4 LiSO4- CaSO4 MnSO4 FeSO4 HSO4- KHSO4 FeSO4+ Fe(SO4)2- H2SO4 S2O8-2 HSO5- NaHSiO3 HSiO3- SiO2 H4(H2SiO4)4-4 H6(H2SiO4)4-2 H2SiO4-2 SrCO3 Sr+2 SrSO4 SrCl+ SrOH+

Although the above ionic strength is high and the choice of database and activity model was possibly not adequate, this should not matter: the reported ionic strength, if calculated using the conventional definition, should accurately match the manually calculated value. This should happen for any database/activity model.

[dlparkhurst]

Phreeqc calculates the ionic strength based on the calculated distribution of species. llnl.dat is an ion-association data base that accounts for activities in part by the formation of ion pairs and ion complexes. Ion pairs and complexes will cause the value to be lower than a value calculated from the analytical data, which effectively assumes all ions are completely dissociated.

Here for example is the calculated distribution of sulfate in the solution that you provided:

Code: [Select]

S(6)          5.000e-01
   SO4-2           1.364e-01   1.324e-02    -0.865    -1.878    -1.013     (0) 
   NaSO4-          1.192e-01   9.140e-02    -0.924    -1.039    -0.115     (0) 
   KSO4-           7.853e-02   6.022e-02    -1.105    -1.220    -0.115     (0) 
   SrSO4           5.229e-02   5.229e-02    -1.282    -1.282     0.000     (0) 
   MgSO4           5.212e-02   5.212e-02    -1.283    -1.283     0.000     (0) 
   LiSO4-          4.952e-02   3.797e-02    -1.305    -1.421    -0.115     (0) 
   CaSO4           7.915e-03   7.915e-03    -2.102    -2.102     0.000     (0) 
   MnSO4           4.003e-03   4.003e-03    -2.398    -2.398     0.000     (0) 
   FeSO4           8.796e-09   8.796e-09    -8.056    -8.056     0.000     (0) 
   HSO4-           1.225e-09   9.391e-10    -8.912    -9.027    -0.115     (0) 
   KHSO4           4.718e-11   4.718e-11   -10.326   -10.326     0.000     (0) 
   FeSO4+          7.082e-17   5.430e-17   -16.150   -16.265    -0.115     (0) 
   Fe(SO4)2-       3.772e-18   2.892e-18   -17.423   -17.539    -0.115     (0) 
   H2SO4           8.137e-23   8.137e-23   -22.090   -22.090     0.000     (0)
Only about a quarter of the sulfate is calculated to be the bare ion SO4-2. The ionic strength will vary by database depending on the log Ks for the aqueous species defined for each database.

pitzer.dat will give a value closest to what you expected because it has very few ion pairs and complexes; although the strongest pairing, like HCO3- (H+ and CO3-2 pair), are accounted for.  Most of the activity correction in the Pitzer model is in the activity coefficients, which are calculated by numerous specific interaction parameters.

[klk]

David, many apologies. The species below were duplicated in the list of species in the -molalities option, and thus my calculated ionic strength had duplicated m*z^2 terms! I just learned that some species are shown more than one time in the output table for the solution composition (Description of solution).

Code: [Select]

NaB(OH)4
BaB(OH)4+
MgB(OH)4+
CaB(OH)4+
BaCO3
BaCl+
NaBr
KBr
SrCO3
MnCO3
CaCO3
MgCO3
NaHCO3
NaCO3-
MgHCO3+
CaHCO3+
MnHCO3+
FeCO3
FeHCO3+
FeCO3+
CaSO4
CaCl+
CaCl2
NaCl
KCl
SrCl+
LiCl
MgCl+
MnCl+
MnCl3-
FeCl+
FeCl4-2
FeCl2
FeCl2+
FeCl+2
FeCl4-
FeSO4
FeSO4+
Fe(SO4)2-
KSO4-
KHSO4
LiSO4-
MgSO4
MnSO4
NaSO4-
NaHSiO3
SrSO4

After removing them, the reported value and my calculation matches perfectly.

Again, apologies for taking the time to clarify anyway and many thanks.

[dlparkhurst]

Here is a script that calculates the ionic strength:

Code: [Select]

SOLUTION
    temperature 60.0
    pressure 100.0
    units mol/kgw
    pH 7.0 charge
    B  0.10
    Ba 1.20
    Br 0.20
    C  4.00
    Ca 0.50
    Cl 2.00
    Fe 0.70
    K  0.80
    Li 0.60
    Mg 0.40
    Mn 0.60
    Na 2.00
    S  0.50
    Si 0.10
    Sr 1.40
USER_PRINT
10 s = SYS("aq", count , name$ , type$ , moles)
20 for i = 1 to count
30   s$ = SPECIES_FORMULA(name$(i), countf, elt$, coef)
40   z = 0
50   for j = 1 to countf
60     if(elt$(j)="charge") then z = coef(j)
70   next j
80   ionic_strength = ionic_strength + 0.5*z^2*moles(i)
90 next i
100 PRINT "Ionic strength: ", ionic_strength, mu
END