zeolite dissolution

[Akhilp123]

I am trying to simulate the dissolution of a zeolite (Heulandite) and the precipitation of Calcite and Kaolinite using the llnl.dat database. However, the results show no significant changes over time?concentrations and mineral amounts remain steady. It appears that the expected geochemical reactions are not occurring.

I have used both EQUILIBRIUM_PHASES and KINETICS blocks in different attempts, but the issue persists.

To assist in identifying the problem, I have attached my input file.
Your support would be greatly appreciated.

Code: [Select]

SOLUTION_SPECIES
CO3-2 + 2H+ = CO2 + H2O
    log_k     16.681
    delta_h   -5.738 kcal
    -analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 0
    -dw       1.92e-09
    -millero  7.29 0.92 2.07 -1.23 -1.6 0
2CO2 = (CO2)2
    log_k     -1.8
    -analytical_expression 8.68 -0.0103 -2190 0 0 0
    -millero  14.58 1.84 4.14 -2.46 -3.2 0
#SiO2 =  SiO2
#    log_k     0

RATES
calcite
-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 scaling factor
10 rem acid solution parameters
11 a1=0
12 E1=0
13 n1=0
20 rem neutral solution parameters
21 a2=6.59E+04
22 E2=66000
30 rem co2 solution parameters
31 a3=1.04E+09
32 E3=67000
33 n2=1.6
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("calcite")
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  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("HCO3-")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
110 rem do not dissolve more minerals than present
#115 if (moles>M) then moles=M
200 save moles
-end

   Kaolinite
-start
 10 REM PARM(1) = SA (BET surface area/kg water)
 20 REM PARM(2) = k25a (rate constant at 25?C in mol/m?s acid mechanism)
 30 REM PARM(3) = k25n (rate constant at 25?C in mol/m?s neutral mechanism)
 40 REM PARM(4) = k25b (rate constant at 25?C in mol/m?s base mechanism)
 50 REM PARM(5) = Eaa (activation energy in J/mol acid mechanism)
 60 REM PARM(6) = Ean (activation energy in J/mol neutral mechanism)
 70 REM PARM(7) = Eab (activation energy in J/mol base mechanism)
 80 REM PARM(8) = na (order of H+ catalysis acid mechanism)
 90 REM PARM(9) = nb (order of H+ catalysis base mechanism)
100 REM PARM(10) = upper pH limit acid mechanism
110 REM PARM(11) = upper pH limit neutral mechanism
120 IF (SI("Kaolinite") < 0) THEN GOTO 130 ELSE GOTO 310
130 REM Dissolution Mechanism
140 IF (M <= 0) THEN GOTO 310
150 t = 1
160 IF (M0 > 0) THEN t = M / M0
170 IF ((-LA("H+")) < PARM(10)) THEN GOTO 230
180 IF ((-LA("H+")) > PARM(11)) THEN GOTO 270
190 REM neutral mechanism
200 k1 = PARM(3)*EXP((PARM(6) / -8.314) * ((1 / TK) - (1 / 298.15)))
210 rate = k1 * PARM(1) * t * (10^(SI("Kaolinite")) - 1)
220 GOTO 300
230 REM acid mechanism
240 k1 = PARM(2)*EXP((PARM(5) / -8.314) * ((1 / TK) - (1 / 298.15)))
250 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(8))) * (10^(SI("Kaolinite")) - 1)
260 GOTO 300
270 REM base mechanism
280 k1 = PARM(4)*EXP((PARM(7) / -8.314) * ((1 / TK) - (1 / 298.15)))
290 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(9))) * (10^(SI("Kaolinite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end

PHASES
CO2(g)
    CO2 = CO2
    log_k     -1.468
    delta_h   -4.776 kcal
    -analytical_expression 10.5624 -0.023547 -3972.8 0 587460 1.9194e-05
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
Calcite
     CaCO3 = Ca+2 + CO3-2
    log_k     -8.48
    delta_h   -2.297 kcal
    -analytical_expression -171.9065 -0.077993 2839.319 71.595 0 0
    -Vm       36.9 cm3/mol
Heulandite
#        Ba.065Sr.175Ca.585K.132Na.383Al2.165Si6.835O18:6 +8.6600 H+  =  + 0.0650 Ba++ + 0.1320 K+ + 0.1750 Sr++ + 0.3830 Na+ + 0.5850 Ca++ + 2.1650 Al+++ + 6.8350 SiO2 + 10.3300 H2O
        Ba.065Sr.175Ca.585K.132Na.383Al2.165Si6.835O18:6H2O +8.6600 H+  =  + 0.0650 Ba++ + 0.1320 K+ + 0.1750 Sr++ + 0.3830 Na+ + 0.5850 Ca++ + 2.1650 Al+++ + 6.8350 SiO2 + 10.3300 H2O
        log_k           3.3506
-delta_H -97.2942 kJ/mol # Calculated enthalpy of reaction Heulandite
# Enthalpy of formation: -10594.5 kJ/mol
        -analytic -1.8364e+001 2.7879e-002 2.8426e+004 -1.7427e+001 -3.4723e+006
#       -Range:  0-300
Kaolinite
        Al2Si2O5(OH)4 +6.0000 H+  =  + 2.0000 Al+++ + 2.0000 SiO2 + 5.0000 H2O
        log_k           6.8101
-delta_H -151.779 kJ/mol # Calculated enthalpy of reaction Kaolinite
# Enthalpy of formation: -982.221 kcal/mol
        -analytic 1.6835e+001 -7.8939e-003 7.7636e+003 -1.2190e+001 -3.2354e+005
#       -Range:  0-300

SOLUTION 1
    temp      35
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 35
    CO2(g)    100

EQUILIBRIUM_PHASES
Heulandite  0 50

KINETICS 1
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        84.1
    -m0       84.1
    -parms    1160 4.9e-12 6.61e-14 8.9e-18 65900 22200 17900 0.78 -0.47 6 8
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.133
    -m0       0.133
    -parms    50000 1.5e-06 3.69e-05 162.6931
    -tol      1e-08
  -steps       2592000 in 500 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5

INCREMENTAL_REACTIONS True
KNOBS
    -iterations            3000
    -convergence_tolerance 1e-10
    -tolerance             1e-15
    -step_size             1000
    -pe_step_size          10
USER_GRAPH 1
    -axis_titles            "Time_Hours" "Delta_Heulandite" ""
    -chart_title            "Delta Heulandite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN_Delta("Heulandite")
  -end
    -active                 true
USER_GRAPH 2
    -axis_titles            "Time_Hours" "Ca Concentration (mol/kgw)" ""
    -chart_title            "Calcium Dissolved"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y TOT("Mg")
  -end
    -active                 true
USER_GRAPH 3
    -axis_titles            "Time_Hours" "Moles" ""
    -chart_title            "Heulandite moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN("Heulandite ")
  -end
    -active                 true
END

[dlparkhurst]

Think.

If you specify equilibrium with Heulandite with EQUILIBRIUM_PHASES, and you do not specify Heulandite in KINETICS, what do you think would be the value of KIN("Heulandite") and KIN_DELTA("Heulandite"). Also, what would be the SI of Heulandite with time?

In USER_GRAPH 2, you are plotting the concentration of Mg. Where is the Mg in your system?

If you don't know what is happening in your system for 1 step, why run for 500 steps?

[Akhilp123]

Thank you for your reply and the helpful questions.

Initially, I had included Heulandite in the KINETICS block. However, due to the lack of reliable dissolution rate data in the literature, I later switched to using EQUILIBRIUM_PHASES to examine the equilibrium behavior instead.

During the first simulation step, I observed a clear trend in calcium concentration, indicating some chemical adjustment. However, the system stabilizes quickly, and beyond that, further steps do not show significant changes. I now realize that running for 500 steps was unnecessary, and I have reduced the number of steps in my revised input accordingly to better focus on the early-time system behavior.
I also wanted to ask about the behavior of Kaolinite in my system. I am not observing any dissolution or precipitation involving Kaolinite?its amount remains unchanged throughout the simulation.
I?m attaching the updated and cleaned input file with this response.
Thank you once again for your valuable guidance.

Code: [Select]

SOLUTION_SPECIES
CO3-2 + 2H+ = CO2 + H2O
    log_k     16.681
    delta_h   -5.738 kcal
    -analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 0
    -dw       1.92e-09
    -millero  7.29 0.92 2.07 -1.23 -1.6 0
2CO2 = (CO2)2
    log_k     -1.8
    -analytical_expression 8.68 -0.0103 -2190 0 0 0
    -millero  14.58 1.84 4.14 -2.46 -3.2 0
#SiO2 =  SiO2
#    log_k     0

RATES
calcite
-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 scaling factor
10 rem acid solution parameters
11 a1=0
12 E1=0
13 n1=0
20 rem neutral solution parameters
21 a2=6.59E+04
22 E2=66000
30 rem co2 solution parameters
31 a3=1.04E+09
32 E3=67000
33 n2=1.6
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("calcite")
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  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("HCO3-")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
110 rem do not dissolve more minerals than present
#115 if (moles>M) then moles=M
200 save moles
-end

   Kaolinite
-start
 10 REM PARM(1) = SA (BET surface area/kg water)
 20 REM PARM(2) = k25a (rate constant at 25?C in mol/m?s acid mechanism)
 30 REM PARM(3) = k25n (rate constant at 25?C in mol/m?s neutral mechanism)
 40 REM PARM(4) = k25b (rate constant at 25?C in mol/m?s base mechanism)
 50 REM PARM(5) = Eaa (activation energy in J/mol acid mechanism)
 60 REM PARM(6) = Ean (activation energy in J/mol neutral mechanism)
 70 REM PARM(7) = Eab (activation energy in J/mol base mechanism)
 80 REM PARM(8) = na (order of H+ catalysis acid mechanism)
 90 REM PARM(9) = nb (order of H+ catalysis base mechanism)
100 REM PARM(10) = upper pH limit acid mechanism
110 REM PARM(11) = upper pH limit neutral mechanism
120 IF (SI("Kaolinite") < 0) THEN GOTO 130 ELSE GOTO 310
130 REM Dissolution Mechanism
140 IF (M <= 0) THEN GOTO 310
150 t = 1
160 IF (M0 > 0) THEN t = M / M0
170 IF ((-LA("H+")) < PARM(10)) THEN GOTO 230
180 IF ((-LA("H+")) > PARM(11)) THEN GOTO 270
190 REM neutral mechanism
200 k1 = PARM(3)*EXP((PARM(6) / -8.314) * ((1 / TK) - (1 / 298.15)))
210 rate = k1 * PARM(1) * t * (10^(SI("Kaolinite")) - 1)
220 GOTO 300
230 REM acid mechanism
240 k1 = PARM(2)*EXP((PARM(5) / -8.314) * ((1 / TK) - (1 / 298.15)))
250 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(8))) * (10^(SI("Kaolinite")) - 1)
260 GOTO 300
270 REM base mechanism
280 k1 = PARM(4)*EXP((PARM(7) / -8.314) * ((1 / TK) - (1 / 298.15)))
290 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(9))) * (10^(SI("Kaolinite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end

PHASES
CO2(g)
    CO2 = CO2
    log_k     -1.468
    delta_h   -4.776 kcal
    -analytical_expression 10.5624 -0.023547 -3972.8 0 587460 1.9194e-05
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
Calcite
     CaCO3 = Ca+2 + CO3-2
    log_k     -8.48
    delta_h   -2.297 kcal
    -analytical_expression -171.9065 -0.077993 2839.319 71.595 0 0
    -Vm       36.9 cm3/mol
Heulandite
#        Ba.065Sr.175Ca.585K.132Na.383Al2.165Si6.835O18:6 +8.6600 H+  =  + 0.0650 Ba++ + 0.1320 K+ + 0.1750 Sr++ + 0.3830 Na+ + 0.5850 Ca++ + 2.1650 Al+++ + 6.8350 SiO2 + 10.3300 H2O
        Ba.065Sr.175Ca.585K.132Na.383Al2.165Si6.835O18:6H2O +8.6600 H+  =  + 0.0650 Ba++ + 0.1320 K+ + 0.1750 Sr++ + 0.3830 Na+ + 0.5850 Ca++ + 2.1650 Al+++ + 6.8350 SiO2 + 10.3300 H2O
        log_k           3.3506
-delta_H -97.2942 kJ/mol # Calculated enthalpy of reaction Heulandite
# Enthalpy of formation: -10594.5 kJ/mol
        -analytic -1.8364e+001 2.7879e-002 2.8426e+004 -1.7427e+001 -3.4723e+006
#       -Range:  0-300
Kaolinite
        Al2Si2O5(OH)4 +6.0000 H+  =  + 2.0000 Al+++ + 2.0000 SiO2 + 5.0000 H2O
        log_k           6.8101
-delta_H -151.779 kJ/mol # Calculated enthalpy of reaction Kaolinite
# Enthalpy of formation: -982.221 kcal/mol
        -analytic 1.6835e+001 -7.8939e-003 7.7636e+003 -1.2190e+001 -3.2354e+005
#       -Range:  0-300

SOLUTION 1
    temp      35
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 35
    CO2(g)    100

EQUILIBRIUM_PHASES
Heulandite  0 50

KINETICS 1
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        84.1
    -m0       84.1
    -parms    1160 4.9e-12 6.61e-14 8.9e-18 65900 22200 17900 0.78 -0.47 6 8
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.133
    -m0       0.133
    -parms    50000 1.5e-06 3.69e-05 162.6931
    -tol      1e-08
  -steps       2592000 in 500 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5

INCREMENTAL_REACTIONS True
KNOBS
    -iterations            3000
    -convergence_tolerance 1e-10
    -tolerance             1e-15
    -step_size             1000
    -pe_step_size          10
USER_GRAPH 2
    -axis_titles            "Time_Hours" "Ca Concentration (mol/kgw)" ""
    -chart_title            "Calcium Dissolved"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y TOT("Ca")
  -end
    -active                 true


[dlparkhurst]

You have the following in the kaolinite rate definition, which does not allow precipitation.

Code: [Select]

120 IF (SI("Kaolinite") < 0) THEN GOTO 130 ELSE GOTO 310
I think you should use the sum of the rates acid, neutral, and base.

The database phreeqc_rates.dat has features to calculate Palandri and Kharaka rates (RATE_PARAMETERS_PK and Basic function RATE_PK) and Hermanska rates (RATE_PARAMETERS_HERMANSKA and Basic function RATE_HERMANSKA), but there are no predefined RATES definitions.

The database kinec_v3.dat has capabilities to calculate rates based on Carbfix (Hermanska and others), and there are predefined RATES definitions.

Don't use KNOBS unless you need to.

[Akhilp123]

Thank you for your previous suggestions. I have made the recommended changes.However, precipitation of both Calcite and Kaolinite still does not occur in my model. The SI values exceed zero, but the kinetic reactant amounts remain constant, and no precipitation is recorded in the output.I would greatly appreciate it if you could review my input file or advise what else might be preventing precipitation.

Code: [Select]

SOLUTION_SPECIES
CO3-2 + 2H+ = CO2 + H2O
    log_k     16.681
    delta_h   -5.738 kcal
    -analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 0
    -dw       1.92e-09
    -millero  7.29 0.92 2.07 -1.23 -1.6 0
2CO2 = (CO2)2
    log_k     -1.8
    -analytical_expression 8.68 -0.0103 -2190 0 0 0
    -millero  14.58 1.84 4.14 -2.46 -3.2 0

RATES
    Heulandite-Ca# CaAl2Si7O18:6H2O
-start
1 name$ = "Heulandite-Ca"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 100
5 S = PARM(2)*TOT("water")
##-----------------Dissolution and precipitation options---------------------##
100  if (PARM(4) = 0) then  goto 1000 else goto 110
110 if (PARM(4) = 1) Then GoTo 150 else goto 200 #
150 If (SR (name$) < 1) Then GoTo 5000 else GoTO 1000 # warning no dissolution reaction
200 If (SR (name$) > 1) Then GoTo 5000 else GoTO 1000 # warning no precipitation reaction#
##------------------Kinetic calculation---------------------##
        #Parameters
1000 Aa = 2.48e-2 #mol.m-2.s-1
1001 An = 1.39e-5 #mol.m-2.s-1
1002 Ab = 3.5e-6 #mol.m-2.s-1
1003 Ea = 33700 #J.mol-1
1004 En = 44200 #J.mol-1
1005 Eb = 44200 #J.mol-1
1006 R =  8.314 #J.deg-1.mol-1
1008 na = 0.82
1009 nb = -0.2
1011 Sig = 7
#rate equations
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na )* S
2001 rplusn = An* (exp(-En/ (R * Tk))) * S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk)))*(act("H+")^nb )* S
2009 rplus = rplusa + rplusn +rplusb
3000 rate = rplus * (1 - (SR("Heulandite-Ca")^(1/Sig)))
4000 moles = rate * time
5000 save moles
-end

 Calcite   #CaCO3; M 100.0869 g/mol
-start
1 name$ = "Calcite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 100
5 S = PARM(2)*TOT("water")
##-----------------Dissolution and precipitation options---------------------##
100  if (PARM(4) = 0) then  goto 1000 else goto 110
110 if (PARM(4) = 1) Then GoTo 150 else goto 200 #
150 If (SR (name$) < 1) Then GoTo 5000 else GoTO 1000 # warning no dissolution reaction
200 If (SR (name$) > 1) Then GoTo 5000 else GoTO 1000 # warning no precipitation reaction#
##------------------Kinetic calculation---------------------##
#Parameters
1000 Aa =5.625# mol.m-2.s-1
1001 Ac = 62.5 # mol.m-2.s-1
1002 Ea =16000# J/mol
1003 Eac =48000# J/mol
1004 R = 8.314 #J.deg-1.mol-1
1006 Sig = 1 
1007 na =1
1008 kc =160
1009 act_c = act("HCO3-")+act("CO3-2")
1010 carb_term = 1-(kc*act_c)/(1+kc*act_c)
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na)* S
2001 rplusc = Ac* (exp(-Eac/ (R * Tk)))*carb_term
2002 rplus = rplusa + rplusc
4000 rate = rplus * (1 - SR("Calcite")^(1/Sig))
5000 moles = rate * time
6000 save moles
-end
Kaolinite # Al2Si2O5(OH)4; M 258.16 g/mol
-start
1 name$ = "Kaolinite"
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * GFW(PHASE_FORMULA(name$)) else S = m0 * ((m/m0)^(2/3)) * GFW(PHASE_FORMULA(name$)) * PARM(2)
4 GOTO 100
5 S = PARM(2)*TOT("water")
##-----------------Dissolution and precipitation options---------------------##
100  if (PARM(4) = 0) then  goto 1000 else goto 110
110 if (PARM(4) = 1) Then GoTo 150 else goto 200 #
150 If (SR (name$) < 1) Then GoTo 5000 else GoTO 1000 # warning no dissolution reaction
200 If (SR (name$) > 1) Then GoTo 5000 else GoTO 1000 # warning no precipitation reaction#
##------------------Kinetic calculation---------------------##
        #Parameters
1000 Aa = 2.85 #mol.m-2.s-1
1001 An = 4.15e-3 #mol.m-2.s-1
1002 Ab = 2.40e-11 #mol.m-2.s-1
1003 Ea = 73000  #J.mol-1
1004 En = 67000 #J.mol-1
1005 Eb = 61000 #J.mol-1
1006 R =  8.314 #J.deg-1.mol-1
1008 na = 0.45
1010 nb = -0.76
1011 Sig = 2
#rate equations
2000 rplusa = Aa* (exp(-Ea/ (R * Tk)))*(act("H+")^na )* S
2001 rplusn = An* (exp(-En/ (R * Tk))) * S
2002 rplusb = Ab* (exp(-Eb/ (R * Tk)))* (act("H+")^nb) * S
2009 rplus = rplusa + rplusn + rplusb
3000 rate = rplus * (1 -SR("Kaolinite")^(1/Sig))
4000 moles = rate * time
5000 save moles
-end

SOLUTION 1
    temp      70
    pH        8
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 50
    -volume 1
    -temperature 35
    CO2(g)    100
KINETICS 1
Heulandite-Ca
    -formula  CaAl2Si7O18:6H2O  1
    -m        1
    -m0       1
    -parms    0 50 0 2
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0
    -m0       0
    -parms    0 0.5 0 1
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0
    -m0       0
    -parms    0 40.42 0 1
    -tol      1e-08
-steps       25920 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 5000
-cvode true
-cvode_steps 1000
-cvode_order 5
USER_GRAPH 1
    -axis_titles            "Time_Hours" "Delta_Kaolinite" ""
    -chart_title            "Delta Kaolinite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN_Delta("Kaolinite")
  -end
    -active                 true
USER_GRAPH 2
    -axis_titles            "Time_Hours" "Ca Concentration (mol/kgw)" ""
    -chart_title            "Calcium leaching"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y TOT("Ca")
  -end
    -active                 true
USER_GRAPH 3
    -axis_titles            "Time_Hours" "Moles" ""
    -chart_title            "Heulandite-Ca moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN("Heulandite-Ca")
  -end
    -active                 true
END

[dlparkhurst]

I think your surface area variable s is zero.

Use PRINT statements in the RATES definition to debug your rate expression.

The following provides the gram formula weight if the mineral is defined in PHASES, GFW(PHASE_FORMULA$("Heulandite-Ca"))