DATABASE phreeqc.dat                          # for K-feldspar
SOLUTION 1
 temp 10; pH 4; Cl 0.1 charge
EQUILIBRIUM_PHASES 1
 Gibbsite 0 0
# CO2(g) -2.5

KINETICS 1
 K-feldspar
 -m0 2.16                                     # 10% K-fsp, mol/L pore water
 -parms 2.26 0.162                            # A0 in m2, V in L
 Quartz; -m0 103; -parms 22.6  0.162          # A0 in m2, V in dm3
 -steps 4.7e7 in 40                           # 1.5 years
INCREMENTAL_REACTIONS true

RATES
 Quartz
 -start
 10 A0 = parm(1)
 20 V = parm(2)
 22 dif_T = 1/TK - 1/298
 30 rate = 10^-(13.7 + 4700*dif_T) * (1 - SR("Quartz")) * A0/V * (m/m0)^0.67 \
           * (1 + 1500*tot("K"))
 32 if SR("Quartz") > 1 then rate = rate * 0.1
 40 moles = rate * time
 50 save moles
 -end

RATES
 K-feldspar
 -start
                  # specific rates from Table 8.7 in mol/m2/s
                  # parm(1) = A in m2, parm(2) = V in L (recalc's sp. rate to mol/kgw)
 1 A0 = parm(1);       2 V = parm(2)
                  # find activities of inhibiting ions...
 3 a_Al = act("Al+3"); 4 BC = act("Na+") + act("K+") + act("Mg+2") + act("Ca+2")
                  # temp corrected with the Arrhenius eqn, Table 8.8
                  # the difference in temperature, TK gives solution temp in Kelvin...
 10 dif_T = 1/TK - 1/271
                  # rate by H+...
 20 pk_H = 11.7 + 3500 * dif_T
 22 rate_H = 10^-pk_H * act("H+")^0.5 / ((1 + a_Al / 4e-6)^0.4 * (1 + BC / 5e-4)^0.15)
                  # rate by hydrolysis...
 30 pk_w = 14.5 + 2000 * dif_T
 32 rate_w = 10^-pk_w / ((1 + a_Al / 4e-6)^0.14 * (1 + BC / 5e-4)^0.15)
                  # rate by OH-...
 40 pk_OH = 13.1 + 2500 * dif_T
 42 rate_OH = 10^-pk_OH
                  # rate by CO2...
 50 pk_CO2 = 13.0 + 2000 * dif_T
 52 rate_CO2 = 10^-pk_CO2 * (10^SI("CO2(g)"))^0.6 
                  # Sum the rate contributions...
 60 rate = rate_H + rate_w + rate_OH + rate_CO2
                  # normalize to mol/kgw, correct for m/m0 and the approach to equi...
 70 rate = rate * A0 / V * (m/m0)^0.67 * (1 - SR("K-feldspar"))
                  # integrate...
 80 moles = rate * time
 90 save moles
 -end

USER_GRAPH
 -head year Si
 -axis_titles Years "mmol Si/L"
 -start
 10 graph_x Total_time / (365 * 24 * 3600)
 20 graph_y tot("Si") * 1e3
 -end
END