Processes > Oxidation and reduction equilibria

Basic question for pyrite reaction

<< < (2/2)

dlparkhurst:
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:

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.

Sincerely,

Jeonghwan Hwang

dlparkhurst:
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 (https://pubs.usgs.gov/of/2004/1068/) 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 (https://geochemsim.org/rate-scripts/).

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.

=========================================================
DATABASE C:\PHREEQC\database\PHREEQC.dat
PHASES
Calcite
CaCO3 = CO3-2 + Ca+2
-log_k   -8.48
-delta_h 0
Gypsum
CaSO4:2H2O = Ca+2 + SO4-2 + 2H2O
-log_k -4.85
-delta_h 0
Siderite
FeCO3 = Fe+2 + CO3-2
-log_k   -10.80
-delta_h 0
Dolomite
CaMg(CO3)2 = Ca+2 + Mg+2 + 2 CO3-2
-log_k   -17.09 #17.90
-delta_h 0
FeS(ppt)
FeS + H+ = Fe+2 + HS-
-log_k   -3.92
-delta_h 0
Fe(OH)3(a)
Fe(OH)3 + 3 H+ = Fe+3 + 3 H2O
-log_k   -4.89
-delta_h 0
END
RATES
Pyrite
-start
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
-end
EXCHANGE_MASTER_SPECIES
Z   Z-
EXCHANGE_SPECIES
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
SURFACE_MASTER_SPECIES
Mont_wa   Mont_waOH
Mont_wb   Mont_wbOH
SURFACE_SPECIES
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
sELECTED_OUTPUT
-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
USER_PUNCH
-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")
end
SOLUTION 1
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
END
USE Solution 1
EXCHANGE 1
NaZ   1.968
CaZ2   0.246
KZ   0.055
MgZ2   0.109
SURFACE 1
-sites_units absolute
Mont_waOH   0.146#0.0627
Mont_wbOH   0.146#0.0627
no_edl
EQUILIBRIUM_PHASES 1
Quartz  0   3.0378
#K-feldspar   0.3947
#Pyrite   0.0   0.0213
Gypsum   0   0.1876
KINETICS 1
Pyrite
-tol    1e-8
-m0     0.0213
-m      0.0213
-parms  5     1
-steps   40000 sec in 100
End
=========================================================

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.
Sincerely,

Jeonghwan Hwang

dlparkhurst:
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.