Processes > Oxidation and reduction equilibria

Basic question for pyrite reaction

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

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.

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.

Navigation

[0] Message Index

[*] Previous page

Go to full version