Processes > Oxidation and reduction equilibria

Basic question for pyrite reaction

<< < (2/2)

Yes, the difference in pe is due to having 1e-30 O2(aq) in the rate calculation and 1e-70 O2(aq) in the equilibrium phases calculation. In the rate calculation, numerical precision of floating point numbers becomes a factor. Within numerical precision of the numerical method, the results of the two calculations are the same.

I think in this case, it does not matter too much, but in general, oxidation of pyrite often leads to precipitation of an iron-oxyhydroxide, like goethite. Your final solution is oversaturated, so I think goethite should be added to equilibrium phases for completeness. If more oxygen is available, its precipitation will make a large difference and lead to acidic waters. [Also, I think you made a sign error in the log K for Fe(OH)3(a)).]

Jeonghwan Hwang:
Thank you for reply.

I think that the 2 models were numerically same but the results were different.
I can say that model 1 returned to a reducing environments (pe=-2.2) whereas model 2 remained an oxidizing environment (pe=8).

If I want to make the same results, I have to change the floating point numbers in light of your comment.
Is it possible? or I have to take another methods?

Additionally, thank you for checking my missing points.


Jeonghwan Hwang

Sorry, my bad. I did not look at the results carefully enough. The calculations do differ.

The Williamson and Rimstedt rate only considers dissolution of pyrite in the presence of O2(aq); it is not designed to reach pyrite equilibrium. Your calculation does indeed dissolve pyrite until O2(aq) is quite small (1e-30); however, the absence of O2 does not lead to pyrite equilibrium. In the rate calculation, O2(aq) is small, but Fe(3) is still significant and S(-2) is small. The rate calculation would also need to consider dissolution of pyrite by Fe+3:

--- Code: ---FeS2 + 14Fe+3 + 8H2O -> 15Fe+2 + 2SO4-2 + 16H+

--- End code ---

Basically, you need to convert all of the Fe(3) to Fe(2) before S(-2) can accumulate to produce SI(Pyrite) = 0. So, one way to modify the rate equation is to add another term in the rate equation that gives the rate of pyrite dissolution as a function of Fe(3) concentration. You could also add a term that is dependent only on (1-SR("Pyrite")).

Palandri and Kharaka ( produced a large compilation of rate equations. I think they have a more complete rate equation for pyrite that includes Fe+3. ZHU, Chen and his students have translated Palandri and Kharaka rates for PHREEQC (

Jeonghwan Hwang:
Based on the report you recommended, the kinetic model was changed.
The Palandri and Kharaka model finally showed a reducing environment with a pe of -2.7.
We are going to apply this model.

   CaCO3 = CO3-2 + Ca+2
   -log_k   -8.48
   -delta_h 0
   CaSO4:2H2O = Ca+2 + SO4-2 + 2H2O
   -log_k -4.85   
   -delta_h 0
   FeCO3 = Fe+2 + CO3-2
   -log_k   -10.80
   -delta_h 0
   CaMg(CO3)2 = Ca+2 + Mg+2 + 2 CO3-2
   -log_k   -17.09 #17.90
   -delta_h 0   
   FeS + H+ = Fe+2 + HS-
   -log_k   -3.92
   -delta_h 0   
   Fe(OH)3 + 3 H+ = Fe+3 + 3 H2O
   -log_k   -4.89
   -delta_h 0
1 rem unit should be mol,Liter-1 and second-1
2 rem parm(1) is surface area in the unit of m2/L
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals M0 is the initial moles of minerals
5 rem parm(2) is a scaling factor
10 rem acid solution parameters
11 a1=2.82E+02      
12 E1=56900 # J/mol 0824
13 n1=-0.500
14 n3=0.500
30 rem neutral solution parameters
31 a2=2.64E+05      
32 E2=56900 # J/mol 0824
33 n2=0.500
36 rem rate=0 if no minerals and undersaturated
41 if (M <=0) then goto 200
42 if (SI("Pyrite") >= 0) THEN GOTO 200
43 SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451 #J/mol/K
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1*ACT("Fe+3")^n3 
80 Rate2=a2*EXP(-E2/R/TK)*ACT("O2")^n2               
90 Rate=(Rate1+Rate2)*(1-SR("pyrite"))*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
   Z   Z-   
   Z- = Z-
      log_k   0.0
   Z- + Na+ = NaZ
      log_K 0.0
   Z- + K+ = KZ
      log_k 0.6
   2Z- + Ca+2 = CaZ2
      log_k 0.41
   2Z- + Mg+2 = MgZ2
      log_k 0.34      
   Mont_wa   Mont_waOH
   Mont_wb   Mont_wbOH
   Mont_waOH = Mont_waOH
   log_k 0
   Mont_waOH + H+ = Mont_waOH2+
   log_k 4.5
   Mont_waOH = Mont_waO- + H+
   log_k -7.9   
   Mont_wbOH = Mont_wbOH
   log_k 0
   Mont_wbOH + H+ = Mont_wbOH2+
   log_k 6.0
   Mont_wbOH = Mont_wbO- + H+
   log_k -10.5   
   -reset false
   -file Kinetic.txt
   -solution true
   -time true
   -pH true
   -pe true
   -water true
   -molalities H2 O2 
   -totals Fe(2) Fe(3) S(-2) S(6)
   -equilibrium_phases Quartz Gypsum
        -headings Mont_waO- Mont_waOH Mont_waOH2+ Mont_wbO- Mont_wbOH Mont_wbOH2+ NaZ KZ CaZ2 MgZ2 CO2(g) O2(g) Ntg(g) Pyrite
   10 PUNCH MOL("Mont_waO-") MOL("Mont_waOH") MOL("Mont_waOH2+") MOL("Mont_wbO-") MOL("Mont_wbOH") MOL("Mont_wbOH2+") MOL("NaZ") MOL("KZ") MOL("CaZ2") MOL("MgZ2")
   20 PUNCH SR("CO2(g)") SR("O2(g)") SR("Ntg(g)")
   30 PUNCH KIN("Pyrite")
   units mol/L   
   pH   7.2
   pe   -2.422
   Temp   15
   C   2.20E-03 as HCO3
   Ca   2.33E-02
   Cl   1.53E-01
   Fe   3.31E-05
   K   8.75E-04
   Mg   9.30E-03
   Na   8.88E-02
   S   6.80E-03 as SO4
   Si   1.85E-04
   -water 0.43
Equilibrium_phases 1
        CO2(g)          -3.51 
   O2(g)           -0.68
   Ntg(g)          -0.11
Save Solution 1
USE Solution 1
      NaZ   1.968
      CaZ2   0.246
      KZ   0.055
      MgZ2   0.109   
   -sites_units absolute
   Mont_waOH   0.146#0.0627   
   Mont_wbOH   0.146#0.0627   
   Quartz  0   3.0378
   #K-feldspar   0.3947   
   #Pyrite   0.0   0.0213
   Gypsum   0   0.1876 
   -tol    1e-8
   -m0     0.0213
   -m      0.0213
   -parms  5     1
-steps   40000 sec in 100

I want to ask you something.
Compared to the Williamson model I used previously (in PHREEQC DB), the dissolution rate of the Palandri and Kharaka model is very fast.
For example, when reacting for 40,000 seconds, the Palandri and Kharaka model approaches equilibrium, but the Williamson model barely responds.

If both models represent meaningful results despite of different kinetic rate, then, how do you think about modifying the parameters of the Williamson model?

After several attempts, I get similar modeling results when I modify the Parm(3) (exponential factor of O2) of the Williamson model from 0.5 to 0.18. (the time to reach equilibrium takes longer for the Williamson model than Palandri and Kharaka model).

The significant question is whether there is a standard for best kinetic model. I think that this part can be affected by the researcher's subjectivity.

Thanks a lot for your help.

Jeonghwan Hwang

I can't help you here. You will have to evaluate the rates that are best applicable to your system.

Clearly, if you can rely on equilibrium calculations instead of rate calculations you will have a stronger case. It is difficult to determine appropriate rate expressions and surface areas, especially in natural systems.


[0] Message Index

[*] Previous page

Go to full version