Pitzer act. coefficients for carbonate w/ dissociation constants by Millero 2006

[Piet_]

Hello,

the aim of my calculations is to compute the activity coefficients of the carbonate- and calcium-ions in a seawater-solution using the Pitzer model. Also i would like to use the dissociation constants for carbonic acid by Millero et al (2006) (Millero, F. J., Graham, T. B., Huang, F., Bustos-Serrano, H., & Pierrot, D. (2006). Dissociation constants of carbonic acid in seawater as a function of salinity and temperature. Marine Chemistry, 100(1-2), 80-94.)(https://www.sciencedirect.com/science/article/pii/S0304420305001921?casa_token=Rtn9kFMDEYgAAAAA:J89RBdc-ecEkrI73RPnyGeLISIFMJ4zadJvG7DPXFc1UkD87BMZC3fehZRs9DykMl97JgBqDeA) due to the difference in pK2* to the standard dissociations constants by Plummer & Busenberg (1982). The purpose for calculating the acitivity coefficients is of course to determine the saturation state of calcite. I am aware how to calculate the saturation state using the standard configurations in PHREEQC but i have problems implementing the aforementioned functions.
I am using the following input file:

Code: [Select]

SELECTED_OUTPUT 1
    -file                 C:\Users\Piet\Documents\Studium\Masterarbeit\Phreeq C\selected_output_FH1 DFZ_F1 T1_FH2_pitzer_3.2.sel
    -simulation           false
    -state                false
    -solution             true
    -distance             false
    -time                 false
    -step                 false
    -ionic_strength       true
    -water                true
    -charge_balance       true
    -percent_error        true
    -totals               C
    -molalities           HCO3-  CO3-2  CO2  Ca+2
    -activities           CO2  CO3-2  HCO3-  Ca+2
    -saturation_indices   Calcite  Barite  Gypsum  Aragonite
    -gases                CO2(g)
SOLUTION 1 FH1-16
    temp      10.87
    pH        7.64
    pe        -0.669
    redox     pe
    units     mmol/l
    density   1
    Br        0.429
    C         2.028
    Ca        5.775
    Cl        309.861
    Fe        13.36 uMol/l
    K         5.422
    Mg        26.743
    Mn        15.763 uMol/l
    Na        254.026
    Oxg       0 mg/l
    S(6)      15.882
    Si        36.888 uMol/l
    Sr        53.8689 uMol/l
    -water    1 # kg
SOLUTION_SPECIES
HCO3- = CO3-2 + H+
    log_k     -9.18
    -analytical_expression -91.291 0 5143.692 14.613358 0 0
CO2 + H2O = HCO3- + H+
    log_k     -5.941
    -analytical_expression -126.7361 0 6320.813 19.568224 0 0

PITZER
-MacInnes   true
-use_etheta true
-redox      false


The corresponding output file:

Code: [Select]

------------------
Reading data base.
------------------

SOLUTION_MASTER_SPECIES
SOLUTION_SPECIES
PHASES
PITZER
EXCHANGE_MASTER_SPECIES
EXCHANGE_SPECIES
SURFACE_MASTER_SPECIES
SURFACE_SPECIES
END
------------------------------------
Reading input data for simulation 1.
------------------------------------

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\pitzer.dat
SELECTED_OUTPUT 1
    file                 C:\Users\Piet\Documents\Studium\Masterarbeit\Phreeq C\selected_output_FH1 DFZ_F1 T1_FH2_pitzer_3.2.sel
    simulation           false
    state                false
    solution             true
    distance             false
    time                 false
    step                 false
    ionic_strength       true
    water                true
    charge_balance       true
    percent_error        true
    totals               C
    molalities           HCO3-  CO3-2  CO2  Ca+2
    activities           CO2  CO3-2  HCO3-  Ca+2
    saturation_indices   Calcite  Barite  Gypsum  Aragonite
    gases                CO2(g)
SOLUTION 1 FH1-16
    temp      10.87
    pH        7.64
    pe        -0.669
    redox     pe
    units     mmol/l
    density   1
    Br        0.429
    C         2.028
    Ca        5.775
    Cl        309.861
    Fe        13.36 uMol/l
    K         5.422
    Mg        26.743
    Mn        15.763 uMol/l
    Na        254.026
    Oxg       0 mg/l
    S(6)      15.882
    Si        36.888 uMol/l
    Sr        53.8689 uMol/l
    water    1 # kg
SOLUTION_SPECIES
HCO3- = CO3-2 + H+
    log_k     -9.18
    analytical_expression -91.291 0 5143.692 14.613358 0 0
CO2 + H2O = HCO3- + H+
    log_k     -5.941
    analytical_expression -126.7361 0 6320.813 19.568224 0 0
PITZER
macinnes   true
use_etheta true
redox      false
-------------------------------------------
Beginning of initial solution calculations.
-------------------------------------------

Initial solution 1. FH1-16

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Br                4.376e-04   4.376e-04
C                 2.069e-03   2.069e-03
Ca                5.891e-03   5.891e-03
Cl                3.161e-01   3.161e-01
Fe                1.363e-05   1.363e-05
K                 5.530e-03   5.530e-03
Mg                2.728e-02   2.728e-02
Mn                1.608e-05   1.608e-05
Na                2.591e-01   2.591e-01
S(6)              1.620e-02   1.620e-02
Si                3.763e-05   3.763e-05
Sr                5.495e-05   5.495e-05

----------------------------Description of solution----------------------------

                                       pH  =   7.640   
                                       pe  =  -0.669   
      Specific Conductance (µS/cm,  11°C)  = 22977
                          Density (g/cm³)  =   1.01455
                               Volume (L)  =   1.00535
                        Activity of water  =   0.990
                 Ionic strength (mol/kgw)  =   3.896e-01
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   9.293e-04
                       Total CO2 (mol/kg)  =   2.069e-03
                         Temperature (°C)  =  10.87
                  Electrical balance (eq)  =  -1.868e-02
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -2.75
                               Iterations  =  14
                         Gamma iterations  =   4
                      Osmotic coefficient  =   0.89812
                         Density of water  =   0.99962
                                  Total H  = 1.110126e+02
                                  Total O  = 5.557577e+01

----------------------------Distribution of species----------------------------

                                                    MacInnes  MacInnes
                                MacInnes       Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma    cm³/mol

   OH-             2.286e-07   1.373e-07    -6.641    -6.862    -0.221     -4.02
   H+              2.898e-08   2.291e-08    -7.538    -7.640    -0.102      0.00
   H2O             5.551e+01   9.898e-01     1.744    -0.004     0.000     18.02
Br            4.376e-04
   Br-             4.376e-04   2.971e-04    -3.359    -3.527    -0.168     24.12
C(4)          2.069e-03
   CO2             1.604e-03   1.696e-03    -2.795    -2.770     0.024     33.71
   CO3-2           2.657e-04   3.742e-05    -3.576    -4.427    -0.851     (0) 
   MgCO3           1.987e-04   1.987e-04    -3.702    -3.702     0.000    -17.07
   HCO3-           0.000e+00   0.000e+00   -51.458   -51.609    -0.150     (0) 
Ca            5.891e-03
   Ca+2            5.891e-03   1.535e-03    -2.230    -2.814    -0.584    -17.40
Cl            3.161e-01
   Cl-             3.161e-01   2.104e-01    -0.500    -0.677    -0.177     17.93
Fe            1.363e-05
   Fe+2            1.363e-05   3.527e-06    -4.866    -5.453    -0.587    -22.14
K             5.530e-03
   K+              5.530e-03   3.793e-03    -2.257    -2.421    -0.164      8.92
Mg            2.728e-02
   Mg+2            2.708e-02   7.635e-03    -1.567    -2.117    -0.550    -20.25
   MgCO3           1.987e-04   1.987e-04    -3.702    -3.702     0.000    -17.07
   MgOH+           1.469e-07   1.403e-07    -6.833    -6.853    -0.020     (0) 
Mn            1.608e-05
   Mn+2            1.608e-05   4.150e-06    -4.794    -5.382    -0.588    -17.36
Na            2.591e-01
   Na+             2.591e-01   1.868e-01    -0.587    -0.729    -0.142     -1.56
S(6)          1.620e-02
   SO4-2           1.620e-02   2.177e-03    -1.790    -2.662    -0.872     15.13
   HSO4-           4.719e-09   3.241e-09    -8.326    -8.489    -0.163     39.48
Si            3.763e-05
   H4SiO4          3.737e-05   3.912e-05    -4.428    -4.408     0.020     53.55
   H3SiO4-         2.605e-07   1.467e-07    -6.584    -6.834    -0.249     27.93
   H2SiO4-2        1.574e-12   1.600e-13   -11.803   -12.796    -0.993     (0) 
Sr            5.495e-05
   Sr+2            5.495e-05   1.428e-05    -4.260    -4.845    -0.585    -16.54

------------------------------Saturation indices-------------------------------

  Phase               SI** log IAP   log K(284 K,   1 atm)

  Akermanite      -18.46     29.28   47.75  Ca2MgSi2O7
  Anhydrite        -1.40     -5.48   -4.07  CaSO4
  Anthophyllite   -14.10     56.91   71.01  Mg7Si8O22(OH)2
  Antigorite      -24.60    481.91  506.51  Mg48Si34O85(OH)62
  Aragonite         0.90     -7.24   -8.14  CaCO3
  Arcanite         -5.41     -7.50   -2.09  K2SO4
  Artinite        -53.74    -32.95   20.79  Mg2CO3(OH)2:3H2O
  Bischofite       -8.27     -3.50    4.77  MgCl2:6H2O
  Bloedite         -6.57     -8.92   -2.35  Na2Mg(SO4)2:4H2O
  Brucite          -4.79    -15.84  -11.06  Mg(OH)2
  Burkeite        -13.35    -14.12   -0.77  Na6CO3(SO4)2
  Calcite           1.06     -7.24   -8.30  CaCO3
  Carnallite      -10.94     -6.60    4.35  KMgCl3:6H2O
  Celestite        -0.92     -7.51   -6.58  SrSO4
  Chalcedony       -0.68     -4.40   -3.72  SiO2
  Chrysotile       -3.37     30.67   34.04  Mg3Si2O5(OH)4
  CO2(g)           -1.49     -2.77   -1.28  CO2
  Diopside         -5.31     16.82   22.13  CaMgSi2O6
  Dolomite          2.97    -13.79  -16.75  CaMg(CO3)2
  Enstatite        -3.29      8.76   12.05  MgSiO3
  Epsomite         -2.85     -4.81   -1.96  MgSO4:7H2O
  Forsterite       -7.74     21.92   29.66  Mg2SiO4
  Gaylussite       -3.73    -13.15   -9.42  CaNa2(CO3)2:5H2O
  Glaserite        -9.29    -13.32   -4.02  NaK3(SO4)2
  Glauberite       -4.34     -9.60   -5.26  Na2Ca(SO4)2
  Goergeyite       -4.55    -34.89  -30.33  K2Ca5(SO4)6H2O
  Gypsum           -0.88     -5.49   -4.60  CaSO4:2H2O
  H2O(g)           -1.89     -0.00    1.89  H2O
  Halite           -2.95     -1.41    1.54  NaCl
  Hexahydrite      -3.28     -4.81   -1.53  MgSO4:6H2O
  Huntite        -196.91   -185.04   11.87  CaMg3(CO3)4
  Kainite          -7.70     -7.89   -0.19  KMgClSO4:3H2O
  Kalicinite       -4.55    -14.49   -9.94  KHCO3
  Kieserite        -4.73     -4.78   -0.05  MgSO4:H2O
  Labile_S         -8.05    -13.72   -5.67  Na4Ca(SO4)3:2H2O
  Leonhardite      -3.91     -4.80   -0.89  MgSO4:4H2O
  Leonite          -8.32    -12.30   -3.98  K2Mg(SO4)2:4H2O
  Magnesite         1.24     -6.54   -7.78  MgCO3
  MgCl2_2H2O      -19.27     -3.48   15.79  MgCl2:2H2O
  MgCl2_4H2O      -10.75     -3.49    7.26  MgCl2:4H2O
  Mirabilite       -2.24     -4.16   -1.92  Na2SO4:10H2O
  Misenite        -73.04    -83.84  -10.81  K8H6(SO4)7
  Nahcolite        -2.05    -12.80  -10.74  NaHCO3
  Natron           -5.10     -5.93   -0.82  Na2CO3:10H2O
  Nesquehonite     -1.39     -6.56   -5.17  MgCO3:3H2O
  Pentahydrite     -3.52     -4.80   -1.28  MgSO4:5H2O
  Pirssonite       -3.90    -13.13   -9.23  Na2Ca(CO3)2:2H2O
  Polyhalite       -9.50    -23.24  -13.74  K2MgCa2(SO4)4:2H2O
  Portlandite     -11.35    -16.54   -5.19  Ca(OH)2
  Quartz           -0.20     -4.40   -4.20  SiO2
  Schoenite        -7.98    -12.31   -4.33  K2Mg(SO4)2:6H2O
  Sepiolite        -3.05     13.10   16.15  Mg2Si3O7.5OH:3H2O
  Sepiolite(d)     -5.56     13.10   18.66  Mg2Si3O7.5OH:3H2O
  SiO2(a)          -1.61     -4.40   -2.79  SiO2
  Sylvite          -3.83     -3.10    0.73  KCl
  Syngenite        -6.84    -12.98   -6.15  K2Ca(SO4)2:H2O
  Talc             -1.21     21.88   23.09  Mg3Si4O10(OH)2
  Thenardite       -3.89     -4.12   -0.23  Na2SO4
  Trona            -7.30    -18.69  -11.38  Na3H(CO3)2:2H2O

**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm.
  For ideal gases, phi = 1.

------------------
End of simulation.
------------------

------------------------------------
Reading input data for simulation 2.
------------------------------------

-------------------------------
End of Run after 0.038 Seconds.
-------------------------------


As you can see the molality of HCO3- is way off. Also the molalities of CO2 and CO3-2 don't match with the values i calculated manually.

[dlparkhurst]

It looks like you used the pitzer.dat database. If so, then you are using activity coefficients based on the Pitzer parameters in the database, which would be inconsistent with Millero. I think Millero would use activity coefficients of 1.0 with all of the activity coefficient correction included in the log Ks.

PHREEQC is probably the wrong application to try to use Millero, but I think you would need to write an entirely new database. I may be wrong about what is needed, but I think you would include only ion pairs and complexes used by Millero. You would include his equilibrium constants for the relevant pairs and complexes. In addition you would need to define activity coefficients of all species as 1.0. You can do this as follows, for example for CO3-2:

Code: [Select]

SOLUTION_SPECIES
HCO3- = CO3-2 + H+
    -gamma 1e8 0
    log_k     -9.18
    -analytical_expression -91.291 0 5143.692 14.613358 0 0

Unfortunately, there is no way in PHREEQC to make the log Ks a function of salinity, so you would probably have to work at a single, fixed salinity, or have a different database for each salinity; I don't see how PHREEQC would handle a calculation of varying salinity with the Millero approach.

I think the Pitzer model should be accurate for the Ca-C system without modifications.

[Piet_]

Thank you for the answer!
Right, I'm using pitzer.dat. It is not quite clear for me, why the dissociation constants by Millero should be inconsistent with Millero. I reread the paper and the following statements (by Millero), to the best of my knowledge, should allow to combine Millero with Pitzer:

1) "the stoichiometric constants are given by
   K1* = [H+][HCO3-] / [CO2]
   K2* = [H+][CO3-2] / [HCO3-]
   with the concentrations being the total stoichiometric values."
   Therefore one should be able to calculate the speciation in moles

2) "The measured stoichiometric constants (K1* and K2*) for carbonic acid in seawater are related to the thermodynamic
    constants by:
    K1* = K1° aH20 yCO2 /yH yHCO3
    K2* = K2° yHCO3 / yH yCO3
    a: activity; y: activity coefficient
    The activity coefficients of H+, HCO3-, CO3-2 and activities of CO2 and H2O can be determined from ionic interaction models
    (Pitzer, 1991)."

I think it is possible to implement the log K as a function of salinity by slightly adapting the fitting equation by Millero.
The general fitting equation is:
pKi - pKi° = Ai + Bi/T + Ci lnT
therefore
pKi = [Ai + Bi/T + Ci lnT] + pKi°
with
pKi°1=-126.34048+6320.813/T+19.568224*LN(T)
In this function the parameter pKi° is only a function of temperature while the parameters A/B/C are a function of salinity.
By calculating pki + pKi° manually in excel the value can be added to the first term for pKi°:
pKi=(-126.34048+[Ai + Bi/T + Ci lnT])+6320.813/T+19.568224*LN(T)

This is not very convenient cause one would need to insert a new Value A1 in the SOLUTION SPECIES for every SOLUTION but it should work in my opinion.

I assume that the Pitzer model in the standard configuration is not accurate cause Millero shows that especially pk2* (determined in seawater) shows large differences to values determined in artificial seawater. Since the standard (thermodynamic) dossociation constants in PHREEQC by Plummer & Busenberg 1982 are determined in CaCO3-H2O-CO2 solution i assume that they show the same differences. I know that the thermodynamic dissociation constants are given in activitys, not in molalities, but since the underlying reason for the difference of the constants in SW and ASW seems to be unknown i assume that it's not considered in the aqueous model of Plummer & Busenberg.

[dlparkhurst]

PHREEQC expects that you are entering the pKi° constants, and then uses activities based on species molality and activity coefficients to solve the mass action equations like  K1° = a(H+)*a(HCO3-)/a(CO2)a(H2O), where the activity is the product of the concentration and the activity coefficient y(H+)*m(H+). Maybe I am misunderstanding, but if you enter K1* instead of K1° in the SOLUTION_SPECIES definition then you are including the activity coefficients in the log K and in the species activity.

Again, my assumption is that Millero's model includes all of the activity correction in K1*, so for his model PHREEQC would need to solve the equation K1* = m(H+)*m(HCO3-)/m(CO2)a(H2O), which means you would have to suppress the activity coefficients (y = 1), relative to the equation set that PHREEQC solves (it always includes an activity coefficient). Now, maybe I am not understanding which log Ks you are actually entering in SOLUTION_SPECIES, but I think it is incorrect if they contain any activity coefficient (salinity) corrections, while still using Pitzer activity coefficients. I think any mixed model between stoichiometric equilibrium constants and intrinsic equilibrium constants is fraught at best. Which Ks contain activity coefficients and which do not? And how do you accommodate both kinds of mass-action equations?

Perhaps you can explain my wrong-headedness, but I think you need to choose Millero's model with stoichiometric log Ks and molality or PHREEQC's model with intrinsic constants and species activities.

[Piet_]

You are right, essentially i tried to implement the stoichiometric constants in the SOLUTION_SPECIES. Embarrassingly I wasn't aware of the problem regarding the different mass equations. Now it's also clear why you suggested to correct the activity coefficient to y=1 by adding the Debye-Hückel parameters. Now, what I don't understand is which activity coefficient is used by PHREEQC to solve the mass equation. It seems like the used activity coefficient is based on the Debye-Hückel equation and not the Pitzer model, since the Pitzer model requires the molality of the species i to calculate the activity coefficient of species i (if im right?). Therefore PHREEQC should solve the mass equation to the "correct" molality, since activity and molality are the same if y=1, as you suggested. Consequently PHREEQC should calculate the correct Pitzer coefficient, if we let your doubts regarding the compatibility of Millero aside.

Therefore it's puzzling that the output speciation and activities are still of, after adapting the debye-hückel parameters:
Input:

Code: [Select]

SELECTED_OUTPUT 1
    -file                 C:\Users\Piet\Documents\Studium\Masterarbeit\Phreeq C\selected_output_FH1 DFZ_F1 T1_FH2_pitzer_3.2.sel
    -simulation           false
    -state                false
    -solution             true
    -distance             false
    -time                 false
    -step                 false
    -ionic_strength       true
    -water                true
    -charge_balance       true
    -percent_error        true
    -totals               C
    -molalities           HCO3-  CO3-2  CO2  Ca+2
    -activities           CO2  CO3-2  HCO3-  Ca+2
    -saturation_indices   Calcite  Barite  Gypsum  Aragonite
    -gases                CO2(g)
SOLUTION 1 FH1-16
    temp      10.87
    pH        7.64
    pe        -0.669
    redox     pe
    units     mmol/l
    density   1.01434
    Br        0.429
    C         2.028
    Ca        5.775
    Cl        309.861
    Fe        13.36 uMol/l
    K         5.422
    Mg        26.743
    Mn        15.763 uMol/l
    Na        254.026
    Oxg       0 mg/l
    S(6)      15.882
    Si        36.888 uMol/l
    Sr        53.8689 uMol/l
    -water    1 # kg
SOLUTION_SPECIES
HCO3- = CO3-2 + H+
    log_k     -9.18
    -analytical_expression -91.291 0 5143.692 14.613358 0 0
    -gamma    100000000 0
CO2 + H2O = HCO3- + H+
    log_k     -5.941
    -analytical_expression -126.7361 0 6320.813 19.568224 0 0
    -gamma    100000000 0

PITZER
-MacInnes   true
-use_etheta true
-redox      false


Output:

Code: [Select]

------------------
Reading data base.
------------------

SOLUTION_MASTER_SPECIES
SOLUTION_SPECIES
PHASES
PITZER
EXCHANGE_MASTER_SPECIES
EXCHANGE_SPECIES
SURFACE_MASTER_SPECIES
SURFACE_SPECIES
END
------------------------------------
Reading input data for simulation 1.
------------------------------------

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\pitzer.dat
SELECTED_OUTPUT 1
    file                 C:\Users\Piet\Documents\Studium\Masterarbeit\Phreeq C\selected_output_FH1 DFZ_F1 T1_FH2_pitzer_3.2.sel
    simulation           false
    state                false
    solution             true
    distance             false
    time                 false
    step                 false
    ionic_strength       true
    water                true
    charge_balance       true
    percent_error        true
    totals               C
    molalities           HCO3-  CO3-2  CO2  Ca+2
    activities           CO2  CO3-2  HCO3-  Ca+2
    saturation_indices   Calcite  Barite  Gypsum  Aragonite
    gases                CO2(g)
SOLUTION 1 FH1-16
    temp      10.87
    pH        7.64
    pe        -0.669
    redox     pe
    units     mmol/l
    density   1.01434
    Br        0.429
    C         2.028
    Ca        5.775
    Cl        309.861
    Fe        13.36 uMol/l
    K         5.422
    Mg        26.743
    Mn        15.763 uMol/l
    Na        254.026
    Oxg       0 mg/l
    S(6)      15.882
    Si        36.888 uMol/l
    Sr        53.8689 uMol/l
    water    1 # kg
SOLUTION_SPECIES
HCO3- = CO3-2 + H+
    log_k     -9.18
    analytical_expression -91.291 0 5143.692 14.613358 0 0
    gamma    100000000 0
CO2 + H2O = HCO3- + H+
    log_k     -5.941
    analytical_expression -126.7361 0 6320.813 19.568224 0 0
    gamma    100000000 0
PITZER
macinnes   true
use_etheta true
redox      false
-------------------------------------------
Beginning of initial solution calculations.
-------------------------------------------

Initial solution 1. FH1-16

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

Br                4.313e-04   4.313e-04
C                 2.039e-03   2.039e-03
Ca                5.806e-03   5.806e-03
Cl                3.115e-01   3.115e-01
Fe                1.343e-05   1.343e-05
K                 5.451e-03   5.451e-03
Mg                2.688e-02   2.688e-02
Mn                1.585e-05   1.585e-05
Na                2.554e-01   2.554e-01
S(6)              1.597e-02   1.597e-02
Si                3.708e-05   3.708e-05
Sr                5.415e-05   5.415e-05

----------------------------Description of solution----------------------------

                                       pH  =   7.640   
                                       pe  =  -0.669   
      Specific Conductance (µS/cm,  11°C)  = 22684
                          Density (g/cm³)  =   1.01434
                               Volume (L)  =   1.00527
                        Activity of water  =   0.990
                 Ionic strength (mol/kgw)  =   3.840e-01
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   9.074e-04
                       Total CO2 (mol/kg)  =   2.039e-03
                         Temperature (°C)  =  10.87
                  Electrical balance (eq)  =  -1.840e-02
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -2.75
                               Iterations  =  14
                         Gamma iterations  =   4
                      Osmotic coefficient  =   0.89818
                         Density of water  =   0.99962
                                  Total H  = 1.110126e+02
                                  Total O  = 5.557476e+01

----------------------------Distribution of species----------------------------

                                                    MacInnes  MacInnes
                                MacInnes       Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma    cm³/mol

   OH-             2.280e-07   1.373e-07    -6.642    -6.862    -0.220     -4.03
   H+              2.900e-08   2.291e-08    -7.538    -7.640    -0.102      0.00
   H2O             5.551e+01   9.900e-01     1.744    -0.004     0.000     18.02
Br            4.313e-04
   Br-             4.313e-04   2.932e-04    -3.365    -3.533    -0.168     24.12
C(4)          2.039e-03
   CO2             1.585e-03   1.675e-03    -2.800    -2.776     0.024     33.71
   CO3-2           2.596e-04   3.695e-05    -3.586    -4.432    -0.847     (0) 
   MgCO3           1.938e-04   1.938e-04    -3.713    -3.713     0.000    -17.07
   HCO3-           0.000e+00   0.000e+00   -51.465   -51.614    -0.149     (0) 
Ca            5.806e-03
   Ca+2            5.806e-03   1.516e-03    -2.236    -2.819    -0.583    -17.40
Cl            3.115e-01
   Cl-             3.115e-01   2.077e-01    -0.507    -0.683    -0.176     17.93
Fe            1.343e-05
   Fe+2            1.343e-05   3.486e-06    -4.872    -5.458    -0.586    -22.15
K             5.451e-03
   K+              5.451e-03   3.744e-03    -2.264    -2.427    -0.163      8.92
Mg            2.688e-02
   Mg+2            2.669e-02   7.541e-03    -1.574    -2.123    -0.549    -20.26
   MgCO3           1.938e-04   1.938e-04    -3.713    -3.713     0.000    -17.07
   MgOH+           1.453e-07   1.386e-07    -6.838    -6.858    -0.021     (0) 
Mn            1.585e-05
   Mn+2            1.585e-05   4.103e-06    -4.800    -5.387    -0.587    -17.38
Na            2.554e-01
   Na+             2.554e-01   1.843e-01    -0.593    -0.735    -0.142     -1.57
S(6)          1.597e-02
   SO4-2           1.597e-02   2.168e-03    -1.797    -2.664    -0.867     15.11
   HSO4-           4.691e-09   3.226e-09    -8.329    -8.491    -0.163     39.47
Si            3.708e-05
   H4SiO4          3.683e-05   3.853e-05    -4.434    -4.414     0.020     53.55
   H3SiO4-         2.557e-07   1.445e-07    -6.592    -6.840    -0.248     27.92
   H2SiO4-2        1.530e-12   1.576e-13   -11.815   -12.802    -0.987     (0) 
Sr            5.415e-05
   Sr+2            5.415e-05   1.412e-05    -4.266    -4.850    -0.584    -16.55

------------------------------Saturation indices-------------------------------

  Phase               SI** log IAP   log K(284 K,   1 atm)

  Akermanite      -18.49     29.26   47.75  Ca2MgSi2O7
  Anhydrite        -1.41     -5.48   -4.07  CaSO4
  Anthophyllite   -14.19     56.82   71.01  Mg7Si8O22(OH)2
  Antigorite      -25.08    481.43  506.51  Mg48Si34O85(OH)62
  Aragonite         0.89     -7.25   -8.14  CaCO3
  Arcanite         -5.42     -7.52   -2.09  K2SO4
  Artinite        -53.75    -32.96   20.79  Mg2CO3(OH)2:3H2O
  Bischofite       -8.28     -3.51    4.77  MgCl2:6H2O
  Bloedite         -6.59     -8.94   -2.35  Na2Mg(SO4)2:4H2O
  Brucite          -4.79    -15.85  -11.06  Mg(OH)2
  Burkeite        -13.40    -14.17   -0.77  Na6CO3(SO4)2
  Calcite           1.05     -7.25   -8.30  CaCO3
  Carnallite      -10.97     -6.62    4.35  KMgCl3:6H2O
  Celestite        -0.93     -7.51   -6.58  SrSO4
  Chalcedony       -0.68     -4.41   -3.72  SiO2
  Chrysotile       -3.40     30.64   34.04  Mg3Si2O5(OH)4
  CO2(g)           -1.49     -2.78   -1.28  CO2
  Diopside         -5.33     16.80   22.13  CaMgSi2O6
  Dolomite          2.94    -13.81  -16.75  CaMg(CO3)2
  Enstatite        -3.31      8.75   12.05  MgSiO3
  Epsomite         -2.86     -4.82   -1.96  MgSO4:7H2O
  Forsterite       -7.75     21.90   29.66  Mg2SiO4
  Gaylussite       -3.75    -13.18   -9.42  CaNa2(CO3)2:5H2O
  Glaserite        -9.32    -13.34   -4.02  NaK3(SO4)2
  Glauberite       -4.36     -9.62   -5.26  Na2Ca(SO4)2
  Goergeyite       -4.60    -34.94  -30.33  K2Ca5(SO4)6H2O
  Gypsum           -0.89     -5.49   -4.60  CaSO4:2H2O
  H2O(g)           -1.89     -0.00    1.89  H2O
  Halite           -2.96     -1.42    1.54  NaCl
  Hexahydrite      -3.28     -4.81   -1.53  MgSO4:6H2O
  Huntite        -196.95   -185.08   11.87  CaMg3(CO3)4
  Kainite          -7.72     -7.91   -0.19  KMgClSO4:3H2O
  Kalicinite       -4.56    -14.50   -9.94  KHCO3
  Kieserite        -4.74     -4.79   -0.05  MgSO4:H2O
  Labile_S         -8.09    -13.76   -5.67  Na4Ca(SO4)3:2H2O
  Leonhardite      -3.92     -4.80   -0.89  MgSO4:4H2O
  Leonite          -8.34    -12.32   -3.98  K2Mg(SO4)2:4H2O
  Magnesite         1.23     -6.55   -7.78  MgCO3
  MgCl2_2H2O      -19.29     -3.50   15.79  MgCl2:2H2O
  MgCl2_4H2O      -10.77     -3.51    7.26  MgCl2:4H2O
  Mirabilite       -2.26     -4.18   -1.92  Na2SO4:10H2O
  Misenite        -73.10    -83.90  -10.81  K8H6(SO4)7
  Nahcolite        -2.06    -12.81  -10.74  NaHCO3
  Natron           -5.12     -5.95   -0.82  Na2CO3:10H2O
  Nesquehonite     -1.40     -6.57   -5.17  MgCO3:3H2O
  Pentahydrite     -3.52     -4.81   -1.28  MgSO4:5H2O
  Pirssonite       -3.93    -13.16   -9.23  Na2Ca(CO3)2:2H2O
  Polyhalite       -9.54    -23.28  -13.74  K2MgCa2(SO4)4:2H2O
  Portlandite     -11.35    -16.54   -5.19  Ca(OH)2
  Quartz           -0.21     -4.41   -4.20  SiO2
  Schoenite        -8.00    -12.33   -4.33  K2Mg(SO4)2:6H2O
  Sepiolite        -3.08     13.07   16.15  Mg2Si3O7.5OH:3H2O
  Sepiolite(d)     -5.59     13.07   18.66  Mg2Si3O7.5OH:3H2O
  SiO2(a)          -1.61     -4.41   -2.79  SiO2
  Sylvite          -3.84     -3.11    0.73  KCl
  Syngenite        -6.86    -13.01   -6.15  K2Ca(SO4)2:H2O
  Talc             -1.26     21.83   23.09  Mg3Si4O10(OH)2
  Thenardite       -3.90     -4.13   -0.23  Na2SO4
  Trona            -7.33    -18.72  -11.38  Na3H(CO3)2:2H2O

**For a gas, SI = log10(fugacity). Fugacity = pressure * phi / 1 atm.
  For ideal gases, phi = 1.

------------------
End of simulation.
------------------

------------------------------------
Reading input data for simulation 2.
------------------------------------

-------------------------------
End of Run after 0.044 Seconds.
-------------------------------


I tried to change the sign of the log K and the terms in the analytical expression which yields a better molality for HCO3- but wrong molalities for CO3-2 and CO2 however.

[dlparkhurst]

Sorry, I was misleading and did not think it all the way through. You cannot use the pitzer.dat database to set activity coefficients to 1 in the way I suggested. Even with if you set all relevant interaction parameters to zero in the PITZER data blocks, each ion will have an activity coefficient based on the D-H term of the Pitzer formulation.

You could use the method I suggested to set activity coefficients to 1.0 in any of the ion-association databases (phreeqc.dat, llnl.dat, minteq.v4.dat, and others).

But, in short, I don't think it is possible to use PHREEQC to create a hybrid model of Pitzer and Millero without  actually changing the code. Even if you wanted to, I don't see how a hybrid model would work, you cannot change just the inorganic carbon species without affecting the other ions; activity of Ca+2, for instance, relies on Ca-HCO3 interaction coefficients, which need to be zero for HCO3, but would need to be non-zero for Ca+2. I think you would run into similar problems even with an ion-association model. The Ca-HCO3- ion pair would probably be removed for Millero, but would then affect the calculated activity of Ca+2. So, I think the best you could do with PHREEQC is to create a full Millero approach by writing a new ion association model database with unit activity coefficients and stoichiometric constants, but it would still only apply at a single salinity. Not very satisfactory, and a fair amount of work.

Note that Millero's model is really only applicable to seawater, whereas Pitzer and the ion-association approach are intended for a wide range of media.

[Piet_]

Okay, bummer. I underestimated the mechanisms one has to consider adding new data in PHREEQC. Seems like my approach wasn't really thought-out.
My samples are a mixture of seawater and freshwater, and since Millero states the equations should be valid for most estuarine and marine waters, i assume they are fairly accurate. I calculated the saturation state for calcite using phreeqc.dat and pitzer.dat in the standard configurations and also using Millero 2006 and Mucci 1983 (Mucci, A. (1983). The solubility of calcite and aragonite in seawater at various salinities, temperatures, and one atmosphere total pressure. Am. J. Sci, 283(7), 780-799.). The Problem is that the total ion activity coefficient products (yT(ca2+)yT(CO3-2) estimated by each method show large differences. The total ion activity coefficient product I calculated for the stoichiometric approach by using the following term:

(yT(ca2+)yT(CO3-2)=(K°/K*)*θ^2

is in the range of 0.021 to 0.017 which agrees with the values calculated by Mucci (1983). The Values calculated by PHREEQC are higher by a factor 2 for pitzer.dat and 4 for phreeqc.dat. To some degree this might be expected due to the lower Mg/Ca ratio in my samples compared to the initial solution Mucci (1983) used to determine the solubility product. But most important, the saturations calculated by PHREEQC, both pitzer.dat and phreeqc.dat, are roughly 60 % less than the saturations I calculated using the stoichiometric approach due to the lower CO3-2 concentrations. This problem was the initial motive to try to implement Milleros equations into PHREEQC.

[dlparkhurst]

I guess it is possible the two approaches differ by that much, but I would have guessed the results would be closer. I wonder if there are some differences in the handling of units, alkalinity, molality of HCO3/CO3, equilibrium constant for calcite, or activity-coefficients calculations (total vs individual ion).

If you post the analytical data that you used and the PHREEQC input file, I will look at the PHREEQC side of the calculations. The saturation index of calcite at 25C in seawater should be the easiest to compare. I would not think it differs by a factor of  2.

[Piet_]

Thank you for the offer. Below you can find the two PHREEQC files which include all the analytical data for one of my samples:

phreeqc.dat:

Code: [Select]

SELECTED_OUTPUT 1
    -file                 C:\Users\Piet\Documents\Studium\Masterarbeit\Phreeq C\selected_output_FH1 DFZ_FH1_16.sel
    -simulation           false
    -state                false
    -solution             true
    -distance             false
    -time                 false
    -step                 false
    -ionic_strength       true
    -water                true
    -charge_balance       true
    -percent_error        true
    -molalities           S-2  HS-  H2S  HCO3-
                          CO3-2  CO2  Ca+2
    -activities           CO2  CO3-2  H2S  HCO3-
                          HS-  S-2  Ca+2
    -saturation_indices   Calcite  Barite  Gypsum  Siderite
                          Rhodochrosite  Sulfur  Aragonite  Pyrite
                          Goethite  Fe(OH)3(a)
    -gases                CO2(g)
    -kinetic_reactants    Calcite
SOLUTION_SPREAD
    -units    mmol/l
 Description O(0)      pe     pH Temperature    N(3)   N(-3)   S(-2)     Ca      K      Mg       Na      Ba      Fe      Mn       P      Si      Sr      F       Cl    S(6)     Br    N(5)      C
            mg/l                           uMol/l uMol/l uMol/l                             uMol/l uMol/l uMol/l uMol/l uMol/l uMol/l                             uMol/l      
      FH1-16     0 -0.669 7.640        10.87   0.336 11.472   0.661 5.775 5.422 26.743 260.116   0.067 13.360 15.763   0.981 36.888 53.869 0.068 309.861 15.882 0.429 21.096 2.028

pitzer.dat:

Code: [Select]

SELECTED_OUTPUT 1
    -file                 C:\Users\Piet\Documents\Studium\Masterarbeit\Phreeq C\selected_output_FH1 DFZ_FH1_16_pitzer.sel
    -simulation           false
    -state                false
    -solution             true
    -distance             false
    -time                 false
    -step                 false
    -ionic_strength       true
    -water                true
    -charge_balance       true
    -percent_error        true
    -totals               C
    -molalities           HCO3-  CO3-2  CO2  Ca+2
    -activities           CO2  CO3-2  HCO3-  Ca+2
    -saturation_indices   Calcite  Barite  Gypsum  Aragonite
    -gases                CO2(g)
    -kinetic_reactants    Calcite
SOLUTION_SPREAD
    -units    mmol/l
 Description O(0)      pe     pH Temperature    N(3)   N(-3)   S(-2)     Ca      K      Mg       Na      Ba      Fe      Mn       P      Si      Sr      F       Cl    S(6)     Br    N(5)      C
            mg/l                           uMol/l uMol/l uMol/l                             uMol/l uMol/l uMol/l uMol/l uMol/l uMol/l                             uMol/l      
      FH1-16     0 -0.669 7.640        10.87   0.336 11.472   0.661 5.775 5.422 26.743 260.116   0.067 13.360 15.763   0.981 36.888 53.869 0.068 309.861 15.882 0.429 21.096 2.028

PITZER
-MacInnes   true
-use_etheta true
-redox      false

The excel file for the manually calculated stoichiometric saturation of the same sample is also included in the attachment.
[Edit: It doesn't seem to be possible to upload excel files to the attachment. Error: "The attachments upload directory is not writable. Your attachment or avatar cannot be saved."]

[dlparkhurst]

I looked at your files. Just a couple small comments.

(1) Pitzer.dat does not have information for F, N(5), N(3), N(-3), O(0), P, or S(-2). They are simply ignored for the calculation.

(2) You have units of per liter, in which case there is a small density correction in converting to moles/kg water. You can either enter the density or "-density 1 calc" in SOLUTION will iterate to find a density for the solution. The difference in mol/kgw should only be 1 or 2 percent.

(3) The ratio to chloride for your solutes seem close to diluted seawater, so I would think the Millero approach should be appropriate. I get a salinity of 19.5 for your water.

(4) phreeqc.dat and pitzer.dat appear to be pretty close in SI(Calcite). I think the only things that could account for a large difference in SI relative to Millero would be the handling of pH and carbon. phreeqc.dat and pitzer.dat assume the pH is the negative log activity of H+. I'm not sure if Millero uses a concentration or activity of H+. As for carbon, you have entered C for PHREEQC, which represents TDIC including CO2(aq). The other option is to enter alkalinity, which would be total alkalinity. If you include boron or phosphorus in the analysis, the associated dissolved species will make a small contribution to the total alkalinity.

(5) Another possibility is the equilibrium constant for calcite. Note that phreeqc.dat and pitzer.dat use different analytical expressions, resulting in slightly different values for the log K calcite.

Still, I just don't see room for a factor of two (0.3 log10 units) error in SI(Calcite) for this system. I'm pretty confident in phreeqc.dat and pitzer.dat. I don't really know how to do the Milllero calculation, but all you would need is the activities of Ca+2 and CO3-2 plus the equilibrium constant (at the correct temperature). At this point, I'll let you decide how you want to proceed.

[Piet_]

Thank you very much for looking into my files and dealing with my problems.

Millero seems to be using total stoichiometric concentrations for H+ (pH on the seawater scale) which I didn't consider. Instead I used the activity. So this might be the reason for the error. Unfortunately I'm not able to correct this mistake because the pH wasn't measured in relation to an appropriate seawater scale, which would be necessary I think.
The H+ concentration should be higher than the activity I guess, so less carbonate, thus a lower SI.