How to prevent a mineral to precipitate?

[Aroua]

Dear PREEQCuSERS?

I developed a code using a kinetic reaction based on Pitzer's model and i want to valided it using Phreeqc as a benchmark!!
In my code i assumed only 4 ionic species to be present in the solution  which are : N+, Cl-, So42-, Ca2+ and after an evaporation through 10 steps a precipitation of only Halite and gypsum will take place.

 It has to be noted that the chemical composition is based on a sea water with rectification of the electroneutrality of the solution  So  please find below the script:

   DATABASE /usr/local/share/doc/phreeqc/database/pitzer.dat
   SOLUTION 1
       temp      25
       pH        7
       pe        4
       redox     pe
       units     mol/kgw
       density   1
       Ca        0.1614
       Cl        0.5606
       Na        0.4928
       S(6)      0.0277 as SO4
       water    1 # kg
   EQUILIBRIUM_PHASES
    Gypsum 0 0
     Anhydrite -1 0         # Mineral_name, SI, Amount
    Halite     0.0  0.0
REACTION #modelling evaporation by removing water from the sample
   H2O -1
   50 moles in 9 steps #~55 moles H2O in 1 L water
   SELECTED_OUTPUT
   reset false
   file Evaporation Gypsum.txt
   si gypsum halite Anhydrite
   ionic_strength   
   molalities Ca+2 Na+ Cl- SO4-2
   activities Ca+2 Na+ Cl- SO4-2
   END

So my question : how to prevent anhydrite or any other mineral to precipitate in the solution except for halite and gypsum ?
i added this script to the code :
  Anhydrite -1 0         # Mineral_name, SI, Amount
but it did not work!!

Many thanks in advance for any help or suggestion !!

best regards,
Aroua

[dlparkhurst]

If you had set Anhydrite +1 0, it would not precipitate, but more simply, remove it from EQUILIBRIUM_PHASES. Put a "#" at the beginning of the line.

[Aroua]

thank you for your response, i've got the result below.
in the initial solution the Anhydrite saturation index is positive and in each subsequent simulation step is negative.
does it mean that Anhydrite precipitate?
and how to eliminate the effect of HSO4- in the activity coefficient.


-------------------------------------------
Beginning of initial solution calculations.
-------------------------------------------

Initial solution 1.   

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

   Elements           Molality       Moles

   Ca                1.614e-01   1.614e-01
   Cl                5.606e-01   5.606e-01
   Na                4.928e-01   4.928e-01
   S(6)              2.770e-02   2.770e-02

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

                                       pH  =   7.000   
                                       pe  =   4.000   
      Specific Conductance (µS/cm,  25°C)  = 57015
                          Density (g/cm³)  =   1.02877
                               Volume (L)  =   1.01124
                        Activity of water  =   0.980
                 Ionic strength (mol/kgw)  =   9.049e-01
                       Mass of water (kg)  =   1.000e+00
                 Total alkalinity (eq/kg)  =   2.545e-08
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =   1.996e-01
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  13.94
                               Iterations  =   7
                         Gamma iterations  =   4
                      Osmotic coefficient  =   0.90052
                         Density of water  =   0.99704
                                  Total H  = 1.110124e+02
                                  Total O  = 5.561702e+01

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

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

   OH-             2.025e-07   9.920e-08    -6.693    -7.004    -0.310     -2.25
   H+              1.408e-07   1.000e-07    -6.851    -7.000    -0.149      0.00
   H2O             5.551e+01   9.800e-01     1.744    -0.009     0.000     18.07
Ca            1.614e-01
   Ca+2            1.614e-01   2.590e-02    -0.792    -1.587    -0.795    -16.53
Cl            5.606e-01
   Cl-             5.606e-01   4.188e-01    -0.251    -0.378    -0.127     18.90
Na            4.928e-01
   Na+             4.928e-01   3.026e-01    -0.307    -0.519    -0.212     -0.36
S(6)          2.770e-02
   SO4-2           2.770e-02   3.041e-03    -1.558    -2.517    -0.959     19.01
   HSO4-           3.631e-08   2.897e-08    -7.440    -7.538    -0.098     41.07

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

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

  Anhydrite         0.15     -4.10   -4.25  CaSO4
  Glauberite       -2.31     -7.66   -5.35  Na2Ca(SO4)2
  Gypsum            0.48     -4.12   -4.60  CaSO4:2H2O
  H2O(g)           -1.51     -0.01    1.50  H2O
  Halite           -2.48     -0.90    1.58  NaCl
  Labile_S         -5.56    -11.23   -5.67  Na4Ca(SO4)3:2H2O
  Mirabilite       -2.40     -3.64   -1.24  Na2SO4:10H2O
  Portlandite     -10.40    -15.59   -5.19  Ca(OH)2
  Thenardite       -3.25     -3.56   -0.30  Na2SO4

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

     5.556e+00 moles of the following reaction have been added:

                    Relative
   Reactant            moles

   H2O                 -1.00000

                    Relative
   Element             moles
   H                   -2.00000
   O                   -1.00000

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

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

Gypsum            0.00    -4.60     -4.60    0.000e+00   1.966e-02   1.966e-02
Halite           -2.39    -0.81      1.58    0.000e+00           0   0.000e+00

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

   Elements           Molality       Moles

   Ca                1.576e-01   1.417e-01
   Cl                6.234e-01   5.606e-01
   Na                5.480e-01   4.928e-01
   S                 8.940e-03   8.039e-03

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

                                       pH  =   6.971      Charge balance
                                       pe  =   4.000      Adjusted to redox equilibrium
      Specific Conductance (µS/cm,  25°C)  = 60357
                          Density (g/cm³)  =   1.02939
                               Volume (L)  =   0.91011
                        Activity of water  =   0.978
                 Ionic strength (mol/kgw)  =   9.189e-01
                       Mass of water (kg)  =   8.992e-01
                 Total alkalinity (eq/kg)  =   2.830e-08
                         Temperature (°C)  =  25.00
                  Electrical balance (eq)  =   1.996e-01
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =  14.75
                               Iterations  =   8
                         Gamma iterations  =   3
                      Osmotic coefficient  =   0.91330
                         Density of water  =   0.99704
                                  Total H  = 9.982268e+01
                                  Total O  = 4.994350e+01

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

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

   OH-             1.875e-07   9.253e-08    -6.727    -7.034    -0.307     -2.23
   H+              1.466e-07   1.070e-07    -6.834    -6.971    -0.137      0.00
   H2O             5.551e+01   9.782e-01     1.744    -0.010     0.000     18.07
Ca            1.576e-01
   Ca+2            1.576e-01   2.651e-02    -0.802    -1.577    -0.774    -16.52
Cl            6.234e-01
   Cl-             6.234e-01   4.666e-01    -0.205    -0.331    -0.126     18.91
Na            5.480e-01
   Na+             5.480e-01   3.358e-01    -0.261    -0.474    -0.213     -0.35
S(6)          8.940e-03
   SO4-2           8.940e-03   9.889e-04    -2.049    -3.005    -0.956     19.05
   HSO4-           1.263e-08   1.008e-08    -7.899    -7.997    -0.098     41.07

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

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

  Anhydrite        -0.33     -4.58   -4.25  CaSO4
  Glauberite       -3.18     -8.53   -5.35  Na2Ca(SO4)2
  Gypsum            0.00     -4.60   -4.60  CaSO4:2H2O
  H2O(g)           -1.51     -0.01    1.50  H2O
  Halite           -2.39     -0.81    1.58  NaCl
  Labile_S         -6.83    -12.51   -5.67  Na4Ca(SO4)3:2H2O
  Mirabilite       -2.81     -4.05   -1.24  Na2SO4:10H2O
  Portlandite     -10.45    -15.64   -5.19  Ca(OH)2
  Thenardite       -3.65     -3.95   -0.30  Na2SO4

[dlparkhurst]

No, anhydrite did not react; it would need to be included in EQUILIBRIUM_PHASES to react. Gypsum did precipitate and it is the stable phase--anhydrite SI is negative, while Gypsum is in equilibrium.

I'm not going to tell you how to remove HSO4- because I don't think you should. Its concentration is negligible in any case.

[Aroua]

Many thanks for this clarification, the result of my code and phreeqc are almost the same for Na, cl ,ca but
not for so4 when i debug phreeqc code source i found almost the same pitzer activity coefficient  but lately
in the code (in reset function class model) phreeqc uses an array named delta to update pitzer activity
i don't want to include this delta in the activity . is it possible
here is a code
x->master[0]->s->la += d;

[dlparkhurst]

I could not find that exact string, but it looks like the way that PHREEQC sequentially updates the master variables to arrive at a solution to the non-linear equations. You definitely don't want to change reset and that basic part of the numerical method. If the program does not converge, that is one thing (and I will lood at it), but if the program does converge, then the results should be the correct solution to the non-linear equations. If the solution is incorrect, please let me know.

HSO4- is orders of magnitude smaller than SO4-2 in the printout that you showed, so its presence should not affect the result significantly.

If you want to reduce HSO4- even further, you can decrease the log K for the reaction. Including the following in your input file will effectively eliminate the effect of HSO4-:

Code: [Select]

SO4-2 + H+ = HSO4-
#log_k    1.979
        log_k -20
#delta_h 4.91    kcal
#-analytic   -5.3585   0.0183412   557.2461
-dw 1.33e-9
-Vm 8.2 9.2590   2.1108   -3.1618 1.1748  0 -0.3 15 0 1 # ref. 1

If you are trying to reproduce results from your code, I don't think the problem is in HSO4-.  I think it is more likely that the Pitzer coefficients are different, or the calculation of non-symmetric mixing factors or some other calculations are different. To debug that, I would probably start with a pure Na2SO4 solution, setting all the PITZER parameters involving SO4-2 in that solution to zero in a PITZER data block in the input file. I would do the same for your code. You can then add parameters to your code and PHREEQC to see where the calculations diverge.