Basic question for pyrite reaction

[Jeonghwan Hwang]

Hello,
I studied about the O2 consumption by pyrite reation.

I made a simple model as below
The solution 1, 2, and 3 showed same results for pe values.

However, when I used this solution with pyrite, the results of solution 1 was different.
May I ask the reason for this?

One more question, based on following sentence 'If the gas reservoir is finite and the total pressure is constant, as in gas bubbles in estuarine and lake sediments and in groundwater, then a fixed-pressure gas phase should be used.' in EXAMPLE 7.
If I want to simulate the transformation of an unsaturated/oxidized underground environment to a saturated/reduced state, would it be better to apply the fixed_pressure method?

Thank you for reading
Sincerely,

Jeonghwan Hwang

==============================
SOLUTION 1
   units mol/L   
   Temp   25
   pH 7
   pe 10.0 O2(g) 0
end

SOLUTION 2
   units mol/L   
   Temp   25
   pH 7
EQUILIBRIUM_PHASES 1
   O2(g)   0
END      

SOLUTION 3
   units mol/L   
   Temp   25
   pH 7
GAS_PHASE 1
        -fixed_pressure
        -pressure       1
        O2(g)          0.21965
END   

#limited O2 Consumption by pyrite

USE Solution 1
EQUILIBRIUM_PHASES 1
   Pyrite   0   1
END      

USE Solution 2
EQUILIBRIUM_PHASES 1
   Pyrite   0   1
END      

USE Solution 3
EQUILIBRIUM_PHASES 1
   Pyrite   0   1
END      

[dlparkhurst]

I think your code did not do quite what you had  in mind for the second two solutions. When you defined SOLUTION 2 ...; EQUILIBRIUM_PHASES 1 ...; END, you are defining a reaction between solution 2 and the equilibrium phases, which react to equilibrium.

The reacted solution has a partial pressure of O2(g) of 1 atm; however, this solution needs to be SAVEd to have it available to react with pyrite. Otherwise, solution 2 remains pure water with a small O2(aq) concentration.

Similarly, SOLUTION 3 reacts with the GAS_PHASE, the bubble of O2(g) partially dissolves to produce a solution with a partial pressure of 1 atm O2(g) and a smaller gas phase volume, but the solution is not saved, so solution 3 remains pure water with very small O2(aq) concentration.

I have added the SAVE statements below. Now all solutions 1-3 have 1 atm partial pressure of O2(g) before they react with pyrite. They each react with the same amount of pyrite to produce pyrite equilibrium.

Code: [Select]

SOLUTION 1
   units mol/L   
   Temp   25
   pH 7
   pe 10.0 O2(g) 0
end

SOLUTION 2
   units mol/L   
   Temp   25
   pH 7
EQUILIBRIUM_PHASES 1
   O2(g)   0
SAVE solution 2
END     

SOLUTION 3
   units mol/L   
   Temp   25
   pH 7
GAS_PHASE 1
        -fixed_pressure
        -pressure       1
        O2(g)          0.21965
END   

#limited O2 Consumption by pyrite

USE Solution 1
EQUILIBRIUM_PHASES 1
   Pyrite   0   1
END     

USE Solution 2
EQUILIBRIUM_PHASES 1
   Pyrite   0   1
END     

USE Solution 3
EQUILIBRIUM_PHASES 1
   Pyrite   0   1
END     

[Jeonghwan Hwang]

Thank you for reply.

I have one more question.
I want to simulate the equilibrium between groundwater and atmosphere, as below.

I think that the Solution 1 (Use Gas_phases) refers the closed system and Solution 2 (use Equilibrium_phases) refers an open system.
Is it correct?

SOLUTION 1
   units mol/L   
   pH   7.2
   pe   -2.42
   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
GAS_PHASE 1
        -fixed_pressure
        -pressure 1
        CO2(g)          0.000308
   O2(g)           0.2095
   Ntg(g)          0.7808
END

SOLUTION 2
   units mol/L   
   pH   7.2
   pe   -2.42
   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
Equilibrium_phases 1
        CO2(g)          -3.51
   O2(g)           -0.68
   Ntg(g)          -0.11
END


Thank you for advice.
Sincerely,

Jeonghwan Hwang

[dlparkhurst]

System 1 has a fixed volume of gas initially. The gases will react with the solution and the equlibrium system will have a different volume of gas and different ratios of gas components than the initial atmospheric ratios.

System 2 will maintain constant partial pressures of each of the gas components.

[Jeonghwan Hwang]

Thanks for your help.

Thanks to this, I got very close to the model I wanted.

1. I used 'Equilibrium_phases' to create an equilibrium reaction of pyrite, gypsum, quartz, and groundwater in oxidizing environment (attachment file 'Pyrite_Equilibrium').

In this model, groundwater initially reacted with the atmosphere (Ntg, O2, CO2) has a high pe (pe>13, oxidizing environment). Next, after reacting with pyrite, it was confirmed that it returned to pe=-2.2 (reducing environment) and O2 was fully reacted with pyrite.

=================================================================
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      
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 Equilibrium.txt
   -solution true
   -time true
   -pH true
   -pe true
   -water true
   -equilibrium_phases Quartz Pyrite 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)
   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)")
   end   
SOLUTION 1
   units mol/L   
   pH   7.2
   pe   -2.42
   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
   Mont_wbOH   0.146
   no_edl   
EQUILIBRIUM_PHASES 1
   Quartz  0   3.03708
   Pyrite   0.0   0.0213
   Gypsum   0   0.1876 
END         
=================================================================

2. Next, I tried to build this model in kinetic ('Pyrite_kinetic')
Based on the Pyrite equation of PHREEQC DB, the  Parm(1) value was changed from 0.3 to 0.75 (log10(0.047[m2/g]*119.85[g/mol])) and modeling was performed.

=================================================================
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
RATES
Pyrite
  -start
1   REM        Williamson and Rimstidt, 1994
2   REM        PARM(1) = log10(specific area), log10(m^2 per mole pyrite)
3   REM        PARM(2) = exp for (M/M0)
4   REM        PARM(3) = exp for O2
5   REM        PARM(4) = exp for H+
10  REM Dissolution in presence of DO
20  if (M <= 0) THEN GOTO 200
30  if (SI("Pyrite") >= 0) THEN GOTO 200
40  log_rate = -8.19 + PARM(3)*LM("O2") + PARM(4)*LM("H+")
50  log_area = PARM(1) + LOG10(M0) + PARM(2)*LOG10(M/M0)
60  moles = 10^(log_area + log_rate) * TIME * (1-SR("Pyrite"))
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_waOHp
   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
   -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
   Mont_wbOH   0.146
   no_edl   
EQUILIBRIUM_PHASES 1
   Quartz  0   3.0378
   #Pyrite   0.0   0.0213
   Gypsum   0   0.1876 
KINETICS 1
Pyrite
   -tol    1e-8
   -m0     0.0213
   -m      0.0213
   -parms  0.75     0.67     0.5      -0.11
   -time 1000 year in 10
End
=================================================================

I think that the results of 1 and 2 should expected to be the same after sufficient times. However, in model 2, the pe was maintained around 8 even after a sufficient time. Looking at the precise results, at 1, the amount of O2 decreased to 10^-68 mol, but at 2 it was about 10^-30 mol. I think that this is the reason.
I want to ask if this reason is correct.

Also, in model 2, I wrote the 'surface complexes', 'cation exchange', 'equilibrium_phases' and 'kinetics' to run simultaneously to implement a complex system. I'd like to ask if there are any potential problems with this approach.

Thanks for reading
Sincerely,

Jeonghwan Hwang

[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: [Select]

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

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.