Strange mineral precipitation behavior?

[ellisone]

Hello,
I am working on some models of mineral dissolution/precipitation, and I have noticed that when using REACTION blocks, I see some mineral precipitation behavior that I wouldn't expect. In particular, if I initialize the model without the elements needed to precipitate the mineral, use EQUILIBRIUM_PHASES to add a mineral, and then start a REACTION that may also be unrelated to the mineral, I sometimes see net precipitation of the mineral when this shouldn't be possible because the required elements aren't provided in excess to allow precipitation.

Below is a minimal example. Reaction Step 1 looks pretty much as I would expect, but in reaction step 2 the delta for goethite is 9.995e-11, indicating that goethite has net-precipitated from solution, and also in an amount that exceeds the Fe that was in solution after step 1.

This seems like an impossible result to me, because there isn't 1e-10 moles of Fe in solution to support this goethite precipitation. Is it just rounding error? Should I not look for this much precision in the delta value?

Thanks a lot for any help!
Eric

   Input file: K:\Raman\Eric\phreeqci\test.pqi
  Output file: K:\Raman\Eric\phreeqci\test.pqo
Database file: C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\phreeqc.dat

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

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

   DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\phreeqc.dat
   SOLUTION 1
       temp      25
       pH        6 charge
       pe        4 O2(g) -0.7
       redox     pe
       units     mmol/l
       density   1
       Fe(3)     0
       P         0
       water    1 # kg
   EQUILIBRIUM_PHASES 1
       Goethite  0 0.55
   REACTION 1
       NaCl  1
       3 millimoles in 3 steps
-------------------------------------------
Beginning of initial solution calculations.
-------------------------------------------

Initial solution 1.   

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

   Elements           Molality       Moles

   Pure water     

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

                                       pH  =   6.997      Charge balance
                                       pe  =  13.625      Equilibrium with O2(g)
      Specific Conductance (µS/cm,  25°C)  = 0
                          Density (g/cm³)  =   0.99704
                               Volume (L)  =   1.00297
                        Activity of water  =   1.000
                 Ionic strength (mol/kgw)  =   1.006e-07
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   6.830e-17
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -6.830e-17
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =   7
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550673e+01

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

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

   OH-             1.006e-07   1.006e-07    -6.997    -6.997    -0.000     -4.14
   H+              1.006e-07   1.006e-07    -6.997    -6.997    -0.000      0.00
   H2O             5.551e+01   1.000e+00     1.744    -0.000     0.000     18.07
H(0)          0.000e+00
   H2              0.000e+00   0.000e+00   -44.394   -44.394     0.000     28.61
O(0)          5.112e-04
   O2              2.556e-04   2.556e-04    -3.592    -3.592     0.000     30.40

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

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

  H2(g)           -41.29    -44.39   -3.10  H2
  H2O(g)           -1.50     -0.00    1.50  H2O
  O2(g)            -0.70     -3.59   -2.89  O2    Pressure   0.2 atm, phi 1.000

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

-----------------------------------------
Beginning of batch-reaction calculations.
-----------------------------------------

Reaction step 1.

Using solution 1.   
Using pure phase assemblage 1.   
Using reaction 1.   

Reaction 1.   

     1.000e-03 moles of the following reaction have been added:

                    Relative
   Reactant            moles

   NaCl                 1.00000

                    Relative
   Element             moles
   Cl                   1.00000
   Na                   1.00000

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Goethite          0.00    -1.00     -1.00    5.500e-01   5.500e-01  -5.017e-14

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

   Elements           Molality       Moles

   Cl                1.000e-03   1.000e-03
   Fe                5.010e-14   5.010e-14
   Na                1.000e-03   1.000e-03

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

                                       pH  =   6.997      Charge balance
                                       pe  =  13.625      Adjusted to redox equilibrium
      Specific Conductance (µS/cm,  25°C)  = 123
                          Density (g/cm³)  =   0.99709
                               Volume (L)  =   1.00299
                        Activity of water  =   1.000
                 Ionic strength (mol/kgw)  =   1.000e-03
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   5.017e-14
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -4.557e-14
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =   8
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550673e+01

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

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

   OH-             1.042e-07   1.005e-07    -6.982    -6.998    -0.016     -4.11
   H+              1.042e-07   1.007e-07    -6.982    -6.997    -0.015      0.00
   H2O             5.551e+01   1.000e+00     1.744    -0.000     0.000     18.07
Cl            1.000e-03
   Cl-             1.000e-03   9.649e-04    -3.000    -3.016    -0.016     18.07
   FeCl+2          3.427e-24   2.976e-24   -23.465   -23.526    -0.061     (0) 
   FeCl+           3.500e-26   3.378e-26   -25.456   -25.471    -0.015     (0) 
   FeCl2+          1.329e-26   1.282e-26   -25.877   -25.892    -0.015     (0) 
   FeCl3           1.237e-30   1.237e-30   -29.908   -29.907     0.000     (0) 
Fe(2)         2.928e-23
   Fe+2            2.917e-23   2.536e-23   -22.535   -22.596    -0.061    -22.14
   FeOH+           8.250e-26   7.964e-26   -25.084   -25.099    -0.015     (0) 
   FeCl+           3.500e-26   3.378e-26   -25.456   -25.471    -0.015     (0) 
   Fe(OH)2         6.729e-30   6.731e-30   -29.172   -29.172     0.000     (0) 
   Fe(OH)3-        2.573e-33   2.483e-33   -32.590   -32.605    -0.015     (0) 
Fe(3)         5.010e-14
   Fe(OH)3         2.753e-14   2.754e-14   -13.560   -13.560     0.000     (0) 
   Fe(OH)2+        2.230e-14   2.153e-14   -13.652   -13.667    -0.015     (0) 
   Fe(OH)4-        2.584e-16   2.494e-16   -15.588   -15.603    -0.015     (0) 
   FeOH+2          7.540e-18   6.547e-18   -17.123   -17.184    -0.061     (0) 
   Fe+3            1.386e-22   1.021e-22   -21.858   -21.991    -0.133     (0) 
   FeCl+2          3.427e-24   2.976e-24   -23.465   -23.526    -0.061     (0) 
   FeCl2+          1.329e-26   1.282e-26   -25.877   -25.892    -0.015     (0) 
   FeCl3           1.237e-30   1.237e-30   -29.908   -29.907     0.000     (0) 
   Fe2(OH)2+4      2.041e-33   1.154e-33   -32.690   -32.938    -0.248     (0) 
   Fe3(OH)4+5      0.000e+00   0.000e+00   -43.898   -44.285    -0.387     (0) 
H(0)          0.000e+00
   H2              0.000e+00   0.000e+00   -44.394   -44.394     0.000     28.61
Na            1.000e-03
   Na+             1.000e-03   9.652e-04    -3.000    -3.015    -0.015     -1.48
   NaOH            9.699e-21   9.701e-21   -20.013   -20.013     0.000     (0) 
O(0)          5.112e-04
   O2              2.556e-04   2.557e-04    -3.592    -3.592     0.000     30.40

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

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

  Fe(OH)3(a)       -5.89     -1.00    4.89  Fe(OH)3
  Goethite          0.00     -1.00   -1.00  FeOOH
  H2(g)           -41.29    -44.39   -3.10  H2
  H2O(g)           -1.50     -0.00    1.50  H2O
  Halite           -7.60     -6.03    1.57  NaCl
  Hematite          2.01     -2.00   -4.01  Fe2O3
  O2(g)            -0.70     -3.59   -2.89  O2

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

Reaction step 2.

Using solution 1.   
Using pure phase assemblage 1.   
Using reaction 1.   

Reaction 1.   

     2.000e-03 moles of the following reaction have been added:

                    Relative
   Reactant            moles

   NaCl                 1.00000

                    Relative
   Element             moles
   Cl                   1.00000
   Na                   1.00000

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Goethite          0.00    -1.00     -1.00    5.500e-01   5.500e-01   9.995e-11

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

   Elements           Molality       Moles

   Cl                2.000e-03   2.000e-03
   Fe                5.042e-14   5.042e-14
   Na                2.000e-03   2.000e-03

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

                                       pH  =   6.997      Charge balance
                                       pe  =  13.625      Adjusted to redox equilibrium
      Specific Conductance (µS/cm,  25°C)  = 244
                          Density (g/cm³)  =   0.99713
                               Volume (L)  =   1.00301
                        Activity of water  =   1.000
                 Ionic strength (mol/kgw)  =   2.000e-03
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   5.049e-14
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -1.184e-13
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =   8
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550673e+01

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

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

   OH-             1.056e-07   1.004e-07    -6.976    -6.998    -0.022     -4.09
   H+              1.056e-07   1.008e-07    -6.976    -6.997    -0.020      0.00
   H2O             5.551e+01   9.999e-01     1.744    -0.000     0.000     18.07
Cl            2.000e-03
   Cl-             2.000e-03   1.903e-03    -2.699    -2.721    -0.022     18.09
   FeCl+2          7.155e-24   5.883e-24   -23.145   -23.230    -0.085     (0) 
   FeCl+           7.011e-26   6.672e-26   -25.154   -25.176    -0.022     (0) 
   FeCl2+          5.251e-26   5.001e-26   -25.280   -25.301    -0.021     (0) 
   FeCl3           9.511e-30   9.516e-30   -29.022   -29.022     0.000     (0) 
Fe(2)         3.097e-23
   Fe+2            3.081e-23   2.540e-23   -22.511   -22.595    -0.084    -22.09
   FeOH+           8.370e-26   7.970e-26   -25.077   -25.099    -0.021     (0) 
   FeCl+           7.011e-26   6.672e-26   -25.154   -25.176    -0.022     (0) 
   Fe(OH)2         6.727e-30   6.730e-30   -29.172   -29.172     0.000     (0) 
   Fe(OH)3-        2.605e-33   2.481e-33   -32.584   -32.605    -0.021     (0) 
Fe(3)         5.042e-14
   Fe(OH)3         2.753e-14   2.754e-14   -13.560   -13.560     0.000     (0) 
   Fe(OH)2+        2.262e-14   2.155e-14   -13.645   -13.667    -0.021     (0) 
   Fe(OH)4-        2.616e-16   2.492e-16   -15.582   -15.603    -0.021     (0) 
   FeOH+2          7.976e-18   6.558e-18   -17.098   -17.183    -0.085     (0) 
   Fe+3            1.554e-22   1.024e-22   -21.808   -21.990    -0.181     (0) 
   FeCl+2          7.155e-24   5.883e-24   -23.145   -23.230    -0.085     (0) 
   FeCl2+          5.251e-26   5.001e-26   -25.280   -25.301    -0.021     (0) 
   FeCl3           9.511e-30   9.516e-30   -29.022   -29.022     0.000     (0) 
   Fe2(OH)2+4      2.559e-33   1.158e-33   -32.592   -32.936    -0.344     (0) 
   Fe3(OH)4+5      0.000e+00   0.000e+00   -43.745   -44.283    -0.538     (0) 
H(0)          0.000e+00
   H2              0.000e+00   0.000e+00   -44.394   -44.394     0.000     28.61
Na            2.000e-03
   Na+             2.000e-03   1.904e-03    -2.699    -2.720    -0.021     -1.46
   NaOH            1.911e-20   1.912e-20   -19.719   -19.718     0.000     (0) 
O(0)          5.112e-04
   O2              2.556e-04   2.557e-04    -3.592    -3.592     0.000     30.40

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

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

  Fe(OH)3(a)       -5.89     -1.00    4.89  Fe(OH)3
  Goethite          0.00     -1.00   -1.00  FeOOH
  H2(g)           -41.29    -44.39   -3.10  H2
  H2O(g)           -1.50     -0.00    1.50  H2O
  Halite           -7.01     -5.44    1.57  NaCl
  Hematite          2.01     -2.00   -4.01  Fe2O3
  O2(g)            -0.70     -3.59   -2.89  O2

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

Reaction step 3.

Using solution 1.   
Using pure phase assemblage 1.   
Using reaction 1.   

Reaction 1.   

     3.000e-03 moles of the following reaction have been added:

                    Relative
   Reactant            moles

   NaCl                 1.00000

                    Relative
   Element             moles
   Cl                   1.00000
   Na                   1.00000

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Goethite          0.00    -1.00     -1.00    5.500e-01   5.500e-01   9.995e-11

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

   Elements           Molality       Moles

   Cl                3.000e-03   3.000e-03
   Fe                5.066e-14   5.066e-14
   Na                3.000e-03   3.000e-03

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

                                       pH  =   6.996      Charge balance
                                       pe  =  13.626      Adjusted to redox equilibrium
      Specific Conductance (µS/cm,  25°C)  = 363
                          Density (g/cm³)  =   0.99717
                               Volume (L)  =   1.00302
                        Activity of water  =   1.000
                 Ionic strength (mol/kgw)  =   3.000e-03
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   5.073e-14
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =  -2.156e-13
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  -0.00
                               Iterations  =   8
                                  Total H  = 1.110124e+02
                                  Total O  = 5.550673e+01

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

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

   OH-             1.066e-07   1.003e-07    -6.972    -6.999    -0.026     -4.08
   H+              1.066e-07   1.009e-07    -6.972    -6.996    -0.024      0.00
   H2O             5.551e+01   9.999e-01     1.744    -0.000     0.000     18.07
Cl            3.000e-03
   Cl-             3.000e-03   2.825e-03    -2.523    -2.549    -0.026     18.10
   FeCl+2          1.108e-23   8.754e-24   -22.955   -23.058    -0.103     (0) 
   FeCl2+          1.172e-25   1.104e-25   -24.931   -24.957    -0.026     (0) 
   FeCl+           1.053e-25   9.920e-26   -24.977   -25.004    -0.026     (0) 
   FeCl3           3.117e-29   3.120e-29   -28.506   -28.506     0.000     (0) 
Fe(2)         3.228e-23
   Fe+2            3.209e-23   2.544e-23   -22.494   -22.594    -0.101    -22.06
   FeCl+           1.053e-25   9.920e-26   -24.977   -25.004    -0.026     (0) 
   FeOH+           8.461e-26   7.976e-26   -25.073   -25.098    -0.026     (0) 
   Fe(OH)2         6.725e-30   6.730e-30   -29.172   -29.172     0.000     (0) 
   Fe(OH)3-        2.630e-33   2.479e-33   -32.580   -32.606    -0.026     (0) 
Fe(3)         5.066e-14
   Fe(OH)3         2.752e-14   2.754e-14   -13.560   -13.560     0.000     (0) 
   Fe(OH)2+        2.287e-14   2.156e-14   -13.641   -13.666    -0.025     (0) 
   Fe(OH)4-        2.640e-16   2.490e-16   -15.578   -15.604    -0.025     (0) 
   FeOH+2          8.318e-18   6.569e-18   -17.080   -17.183    -0.103     (0) 
   Fe+3            1.689e-22   1.026e-22   -21.772   -21.989    -0.216     (0) 
   FeCl+2          1.108e-23   8.754e-24   -22.955   -23.058    -0.103     (0) 
   FeCl2+          1.172e-25   1.104e-25   -24.931   -24.957    -0.026     (0) 
   FeCl3           3.117e-29   3.120e-29   -28.506   -28.506     0.000     (0) 
   Fe2(OH)2+4      3.030e-33   1.161e-33   -32.519   -32.935    -0.416     (0) 
   Fe3(OH)4+5      0.000e+00   0.000e+00   -43.631   -44.281    -0.651     (0) 
H(0)          0.000e+00
   H2              0.000e+00   0.000e+00   -44.394   -44.394     0.000     28.61
Na            3.000e-03
   Na+             3.000e-03   2.827e-03    -2.523    -2.549    -0.026     -1.44
   NaOH            2.835e-20   2.837e-20   -19.547   -19.547     0.000     (0) 
O(0)          5.112e-04
   O2              2.556e-04   2.558e-04    -3.592    -3.592     0.000     30.40

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

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

  Fe(OH)3(a)       -5.89     -1.00    4.89  Fe(OH)3
  Goethite          0.00     -1.00   -1.00  FeOOH
  H2(g)           -41.29    -44.39   -3.10  H2
  H2O(g)           -1.50     -0.00    1.50  H2O
  Halite           -6.67     -5.10    1.57  NaCl
  Hematite          2.01     -2.00   -4.01  Fe2O3
  O2(g)            -0.70     -3.59   -2.89  O2

**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.079 Seconds.
-------------------------------
 

[dlparkhurst]

Looks like a bug, but not one I'm going to worry about. For your case, if the solution does not have any iron in the system, then PHREEQC dissolves 1e-10 moles of goethite to do the calculation. After the numerical solution is found, it tries to adjust the mole transfer for that 1e-10 moles that were added. It works correctly for the first step, but it looks like the subsequent steps do not account for that 1e-10 moles correctly. I think if the mole transfer were greater than 1e-10, there would not be a problem.

If you (1) add some Fe to the initial solution or (2) use INCREMENTAL_REACTIONS, I think you will get correct results. As it is, the error is on the order of 1e-10 moles in the mole transfer of Goethite, but the solution composition should be identical to the results using (1) or (2).