Mineral Dissolution as a Function of pH

[OmranS]

Hi there,

I aim to simulate the dissolution of barite and pyrite. My goal is to see how long it would take for saturation to be reached in solutions with different pHs. The issue I'm running into is saturation is being reached in less than a second no matter what pH or mineral is inputted.

Here is my input for pyrite

Code: [Select]

# Title Pyrite Dissolution

PHASES
# from https://github.com/HydrogeoIU/PHREEQC-Kinetic-Library
Pyrite
Fe1S2 + 1.0000 H2O = 0.5000 O2    + 2.0000 HS- + 1.0000 Fe+2
-analytic 7.990927e+003 2.426724e+000 -3.481962e+005 -3.112050e+003 1.492998e+007 -8.237808e-004
-Vm 23.9400

SOLUTION
temp 25
units mol/kgw
ph 3


EQUILIBRIUM_PHASES
Pyrite 0 1

KINETICS 1
Pyrite
 -m 1           # moles of solid per kg of water
 -parms 50  1       # total surface area per kg of water (m2/kgw) and the scaling factor
 -steps 1 second in 100 steps        # define time steps
 INCREMENTAL_REACTIONS TRUE


RATES
pyrite
# from Palandri and Kharaka 2004
# experimental condition range T=20-40C, pH=1-4

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
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 correction factor
10 rem acid solution parameters
11 a1=2.82E+02
12 E1=56900
13 n1=-0.500
14 n3=0.500
30 rem neutral solution parameters
31 a3=2.64E+05
32 E3=56900
33 n2=0.500
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("pyrite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1*ACT("Fe+3")^n3  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)*ACT("O2")               #neutral rate expression
90 Rate=(Rate1+Rate2)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end


SELECTED_OUTPUT
     -file pyrite_dissolution.txt
     -reset false
USER_PUNCH
     -headings Time(seconds) SI(Pyrite) pH Tot_S Tot_Fe
     -start
     10 punch TOTAL_TIME SI("Pyrite") LA("H+")*(-1) TOT("S") TOT("Fe")
     -end
END

I would appreciate help to find out why this is happening

[dlparkhurst]

By putting Pyrite in EQUILIBRIUM_PHASES (with positive moles), the solution will always be in equilibrium with pyrite, regardless of the KINETICS reaction.

If you want pyrite to react kinetically, you should remove pyrite from EQUILIBRIUM_PHASES. Then you need to fix up your RATES definition. Your pure water solution has no Fe+3 and very little O2(aq), so rate1 and rate2 would be small. You should probably define some O(0) (dissolved O2) in SOLUTION. Further, you need to define A2 and E2 for the rate2 calculation.

[OmranS]

Thank you very much for your reply!

As per your suggestion, pyrite was removed from the equilibrium phase datablock. It was replaced with Goethite, a phase also created in the dissolution process. Oxygen was specified in the solution data block, yet this created an odd dissolution curve when concentration of iron over time was plotted (this plot has been attached)

My input is now:

Code: [Select]

# Title Pyrite Dissolution

PHASES
# from https://github.com/HydrogeoIU/PHREEQC-Kinetic-Library
Pyrite
Fe1S2 + 1.0000 H2O = 0.5000 O2    + 2.0000 HS- + 1.0000 Fe+2
-analytic 7.990927e+003 2.426724e+000 -3.481962e+005 -3.112050e+003 1.492998e+007 -8.237808e-004
-Vm 23.9400

SOLUTION
temp 25
units mol/kgw
O(0) 1.0 O2(g) -0.7
ph 3


EQUILIBRIUM_PHASES
Goethite

KINETICS 1
Pyrite
 -m 1           # moles of solid per kg of water
 -parms 50  1       # total surface area per kg of water (m2/kgw) and the scaling factor
 -steps 518400 second in 100 steps        # 6 days
 INCREMENTAL_REACTIONS TRUE


RATES
pyrite
# from Palandri and Kharaka 2004
# experimental condition range T=20-40C, pH=1-4

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
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 correction factor
10 rem acid solution parameters
11 a1=2.82E+02
12 E1=56900
13 n1=-0.500
14 n3=0.500
30 rem neutral solution parameters
31 a3=2.64E+05
32 E3=56900
33 n2=0.500
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("pyrite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1*ACT("Fe+3")^n3  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)*ACT("O2")               #neutral rate expression
90 Rate=(Rate1+Rate2)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end


SELECTED_OUTPUT
     -file pyrite_dissolution.txt
     -reset false
USER_PUNCH
     -headings Time(seconds) SI(Pyrite) pH Tot_S Tot_Fe
     -start
     10 punch TOTAL_TIME/86400 SI("Pyrite") LA("H+")*(-1) TOT("S") TOT("Fe")
     -end
END


My confusion comes from the reason for the concentration to remain stagnant until the approximate 1.4 days mark, where it suddenly jumps on the attached plot

edit: the figure will not attach and is showing an error. Looks like viewing and adding attachments hasn't worked for a while as said by

Yes, the attachments haven't been functional for a while, both to view old ones or to add new ones.

If you have a particular one you are interested in, I can check whether it can be recreated and added within a post.

So if I were to describe it, the plot follows the x-axis (y=0) until ~1.44 days is reached, where it suddenly jumps to a logarithmic or square root-like trend

Thank you for your time and help

[dlparkhurst]

The jump in saturation index corresponds to the point where O2(aq) is depleted, at which point the activities of Fe+2 and HS- increase dramatically. SI(Pyrite) is directly related to the sum of the three log activities that are plotted in USER_GRAPH. Note that the SI increases, but in fact, the rates (Rate1) is decreasing.

I again note you have an error in Rate2. A3, E3, and n2 defined in lines 12-14 of the RATES definition are not used in the calculation, and A2 and E2, which are used in line 80 are undefined and default to zero. Therefore Rate2 is always zero.

Code: [Select]

# Title Pyrite Dissolution

PHASES
# from https://github.com/HydrogeoIU/PHREEQC-Kinetic-Library
Pyrite
Fe1S2 + 1.0000 H2O = 0.5000 O2    + 2.0000 HS- + 1.0000 Fe+2
-analytic 7.990927e+003 2.426724e+000 -3.481962e+005 -3.112050e+003 1.492998e+007 -8.237808e-004
-Vm 23.9400
RATES
pyrite
# from Palandri and Kharaka 2004
# experimental condition range T=20-40C, pH=1-4

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
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 correction factor
10 rem acid solution parameters
11 a1=2.82E+02
12 E1=56900
13 n1=-0.500
14 n3=0.500
30 rem neutral solution parameters
31 a3=2.64E+05
32 E3=56900
33 n2=0.500
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("pyrite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1*ACT("Fe+3")^n3  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)*ACT("O2")               #neutral rate expression
90 Rate=(Rate1+Rate2)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles

-end
END
SOLUTION
temp 25
units mol/kgw
O(0) 1.0 O2(g) -0.7
ph 3


EQUILIBRIUM_PHASES
Goethite

KINETICS 1
Pyrite
-m 1               # moles of solid per kg of water
-parms 50  1       # total surface area per kg of water (m2/kgw) and the scaling factor
-steps 518400 second in 100 steps        # 6 days
INCREMENTAL_REACTIONS TRUE

USER_GRAPH 1
    -headings               time SI(Pyrite) O2 HS- Fe+2
    -axis_scale sy_axis     auto auto auto auto auto
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/86400
20 GRAPH_Y SI("Pyrite")
30 GRAPH_Y 0.5*LA("O2"), 2*LA("HS-"), LA("Fe+2")
  -end
    -active                 true

END