RATES
 Pyrite  # rates from data compiled by Williamson and Rimstidt 1994, GCA 58, 5443
 -start
  1 A = 120 * m0                                                                                            # initial surface area in m2 indicates 1 um size crystals
 10 if SI("Pyrite")>0 then goto 100                                                                                                     # step out when supersaturated...
 20 fH = mol("H+")
 30 fFe2 = (1 + tot("Fe(2)") / 1e-6)
 40 if mol("O2") < 1e-6 then goto 80
 50 rO2 = 10^-8.19 * mol("O2")^0.5 * fH^-0.11                                                                                                           # ...rate with oxygen
 60 rO2_Fe3 = 6.3e-4 * tot("Fe(3)")^0.92 * fFe2^-0.43                                                                                      # ...rate with oxygen and Fe3+
 70 goto 90
 80 rem                                                                                                                # rate with Fe3+ without oxygen, and for pH < 3...
 81 rFe3 = 1.9e-6 * tot("Fe(3)")^0.28 * fFe2^-0.52 * fH^-0.3
 90 rate = A * (m/m0)^0.67 * (rO2 + rO2_Fe3 + rFe3) * (1 - SR("Pyrite"))                                                              # ...sum terms, zero at equilibrium
 100 save rate * time
 -end

SOLUTION_SPECIES                                                                                                                          # make N2 the only N species...
 2NO3- + 12H+ + 10e- = N2 + 6H2O; log_k 500 #207.080

SOLUTION 1
 -water 0.0069239
 -temp 20
 pH 7    charge;     pe 14  O2(g)  -1.0878
 Ca 1    Calcite;    C  1   CO2(g) -2.6021
 Fe 1e-3 Goethite 2; N  1.3 N2(g)  -0.0382

EQUILIBRIUM_PHASES 1
 Goethite 2; Calcite 0; Gypsum 0 0

KINETICS 1
  Pyrite; -m0 1.32e-4; -step 0 10*5e5
INCREMENTAL_REACTIONS true

GAS_PHASE 1
 -fixed_pressure  # default 1 atm
 -volume 0.02127; -temp 20.0
 CO2(g) 0.0025; O2(g) 0.0817; N2(g) 0.9157

USER_GRAPH
 -headings time CO2 O2 N2
 -axis_titles "Time / hour" "P_O2 and P_CO2 / atm" "P_N2 / atm"
 -start
 10 graph_x total_time / 3600
 20 graph_y 10^si("CO2(g)"), 10^si("O2(g)"); 30 graph_sy 10^si("N2(g)")
 -end
END