CO2 dissolved in water to react with rock

[ahmalisha]

I try to simulate the interaction  between the CO2 dissolved  in the water with the rock and how that affect the porosity of the rock by using the code attached but i have not noticed any change in the porosity or volume 

Code: [Select]

RATES
    Calcite
-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("Calcite") < 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("Calcite")) - 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("Calcite")) - 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("Calcite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
    Dolomite
-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("Dolomite") < 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("Dolomite")) - 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("Dolomite")) - 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("Dolomite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Quartz
-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("Quartz") < 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("Quartz")) - 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("Quartz")) - 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("Quartz")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
K-Feldspar
-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("K-Feldspar") < 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("K-Feldspar")) - 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("K-Feldspar")) - 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("K-Feldspar")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Siderite
-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("Siderite") < 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("Siderite")) - 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("Siderite")) - 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("Siderite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Pyrite
-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("Pyrite") < 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("Pyrite")) - 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("Pyrite")) - 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("Pyrite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 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
Illite
-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("Illite") < 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("Illite")) - 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("Illite")) - 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("Illite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Albite
-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("Albite") < 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("Albite")) - 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("Albite")) - 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("Albite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Smectite-high-Fe-Mg
-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("Smectite-high-Fe-Mg") < 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("Smectite-high-Fe-Mg")) - 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("Smectite-high-Fe-Mg")) - 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("Smectite-high-Fe-Mg")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Gypsum
-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("Gypsum") < 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("Gypsum")) - 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("Gypsum")) - 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("Gypsum")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
END
TITLE CO2-H2S Rock Fluid Interaction
SOLUTION 1
    temp      140
    pH        5.5  # Slightly acidic to drive carbonate dissolution
    pe        4
    redox     pe
    units     mg/l
    C(4)      1000  # Higher initial CO2 to prompt dissolution
    Ca        54.6
    Cl        18927
    Fe        0.02
    K         111
    Mg        12.1
    Na        8504
    S(6)      1300
    Sr        0.05
    -water    1 # kg
GAS_PHASE 1
    -fixed_volume               # Switch to a fixed volume to allow CO? to dissolve
    -volume 1                   # Volume in liters
    CO2(g)    81.5292           # CO? partial pressure in atm
    H2S(g)    0.126             # H?S partial pressure in atm
    -temperature 140            # Temperature in ?C

EQUILIBRIUM_PHASES 1
    Albite                   0 190
    Gypsum                   0 100
    Calcite                  0 78200
    Dolomite                 0 60900
    Illite                   0 14400
    K-feldspar               0 152
    Kaolinite                0 84100
    Pyrite                   0 60800
    Quartz                   0 648400
    Siderite                 0 46600
    Smectite-high-Fe-Mg      0 54.6
SAVE SOLUTION 2
END
USE SOLUTION 2

GAS_PHASE 1
    -fixed_pressure
    -pressure 81.6552
    -volume 1
    -temperature 140
    CO2(g)    810.5292
    H2S(g)    0.126

KINETICS 1
Calcite
    -formula  CaCO3  1
    -m        78.2
    -m0       78.2
    -parms    67.1 0.501 1.55e-06 0.0003311 14400 23500 35400 1 1 6 8
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        60.9
    -m0       60.9
    -parms    63.6 6.457e-06 2.9512e-08 7.76e-06 36100 52200 34800 1 1 6 8
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        648.4
    -m0       648.4
    -parms    68.6 0 1.023e-14 0 0 87700 0 0 0 6 8
    -tol      1e-08
K-feldspar
    -formula  KAlSi3O8  1
    -m        0.152
    -m0       0.152
    -parms    71.1 8.71e-11 3.89e-13 6.31e-22 51700 38000 94100 0.5 0.82 6 8
    -tol      1e-08
Siderite
    -formula  FeCO3  1
    -m        46.6
    -m0       46.6
    -parms    46.2 0.000158 1.26e-09 0 45000 62800 0 0.9 0 6 8
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        60.8
    -m0       60.8
    -parms    36.3 4.89e-12 6.6e-14 8.91e-18 65900 22200 17900 0.78 0 6 8
    -tol      1e-08
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
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        14.4
    -m0       14.4
    -parms    460 1.95e-13 3.9e-15 3.9e-15 48000 48000 48000 0.22 -0.13 6 8
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.19
    -m0       0.19
    -parms    69.4 6.91e-11 2.75e-13 0 65000 69800 0 0.46 0 6 8
    -tol      1e-08
Smectite-high-Fe-Mg
    -formula  Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12  1
    -m        0.0546
    -m0       0.0546
    -parms    1040 1.047e-11 1.659e-13 3.019e-17 23600 0 58900 0.34 -0.4 6 8
    -tol      1e-08
Gypsum
    -formula  CaSO4:2H2O  1
    -m        0.1
    -m0       0.1
    -parms    61.6 0 0.001622 0 0 0 0 0 0 6 8
    -tol      1e-08
-steps       1e0 1e1 1e2 1e3 1e4 1e5 1e6 1e7 6.3072e7 # Seconds up to exactly 2 year
INCREMENTAL_REACTIONS True
USER_GRAPH 1
    -headings               time Calcite Dolomite Quartz K-feldspar \
Siderite Pyrite Kaolinite Illite \
Albite Smectite-high-Fe-Mg Gypsum
    -axis_titles            "Seconds" "Saturation Index" ""
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y SI("Calcite"), SI("Dolomite"), SI("Quartz"), SI("K-feldspar")
30 GRAPH_Y SI("Siderite"), SI("Pyrite"), SI("Kaolinite"), SI("Illite")
40 GRAPH_Y SI("Albite"), SI("Smectite-high-Fe-Mg"), SI("Gypsum")
  -end
USER_GRAPH 2
    -headings               time Calcite Dolomite Quartz K-feldspar \
Siderite Pyrite Kaolinite Illite \
Albite Smectite-high-Fe-Mg Gypsum
    -axis_titles            "Seconds" "Saturation Index" ""
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y KIN_DELTA("Calcite"), KIN_DELTA("Dolomite"), KIN_DELTA("Quartz"), KIN_DELTA("K-feldspar")
30 GRAPH_Y KIN_DELTA("Siderite"), KIN_DELTA("Pyrite"), KIN_DELTA("Kaolinite"), KIN_DELTA("Illite")
40 GRAPH_Y KIN_DELTA("Albite"), KIN_DELTA("Smectite-high-Fe-Mg"), KIN_DELTA("Gypsum")
  -end
USER_GRAPH 2
    -headings               time Calcite Dolomite Quartz K-feldspar \
Siderite Pyrite Kaolinite Illite \
Albite Smectite-high-Fe-Mg Gypsum
    -axis_titles            "Seconds" "Reaction Rates" ""
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
-start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y KIN_DELTA("Calcite"), KIN_DELTA("Dolomite"), KIN_DELTA("Quartz"), KIN_DELTA("K-feldspar")
30 GRAPH_Y KIN_DELTA("Siderite"), KIN_DELTA("Pyrite"), KIN_DELTA("Kaolinite"), KIN_DELTA("Illite")
40 GRAPH_Y KIN_DELTA("Albite"), KIN_DELTA("Smectite-high-Fe-Mg"), KIN_DELTA("Gypsum")
-end

USER_GRAPH 3
    -headings               "Time (Days)" "Volume Change of Calcite" "Volume Change of Dolomite" "Volume Change of Quartz" \
                            "Volume Change of K-feldspar" "Volume Change of Siderite" "Volume Change of Pyrite" \
                            "Volume Change of Kaolinite" "Volume Change of Illite" "Volume Change of Albite" \
                            "Volume Change of Smectite-high-Fe-Mg" "Volume Change of Gypsum"
    -axis_titles            "Time (Days)" "Phase Volume Changes (cm?/100g)"
    -chart_title            "Phase Volume Changes for Each Mineral Over Time"
    -axis_scale             x_axis log
-start
10 GRAPH_X TOTAL_TIME / 86400
20 v_calcite = KIN("Calcite") * 36.89
30 v_dolomite = KIN("Dolomite") * 64.5
40 v_quartz = KIN("Quartz") * 22.69
50 v_kfeldspar = KIN("K-feldspar") * 108.15
60 v_siderite = KIN("Siderite") * 46.2
70 v_pyrite = KIN("Pyrite") * 36.3
80 v_kaolinite = KIN("Kaolinite") * 99.35
90 v_illite = KIN("Illite") * 102.0
100 v_albite = KIN("Albite") * 100.2
110 v_smectite = KIN("Smectite-high-Fe-Mg") * 1040
120 v_gypsum = KIN("Gypsum") * 58.3
130 GRAPH_Y v_calcite, v_dolomite, v_quartz, v_kfeldspar, v_siderite, v_pyrite, v_kaolinite, v_illite, v_albite, v_smectite, v_gypsum
-end
-active true

[dlparkhurst]

Well, you could investigate with one kinetic reaction and look carefully at the output.

I think your mineral volumes are changing, just not very much.

Code: [Select]

RATES
    Calcite
-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("Calcite") < 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("Calcite")) - 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("Calcite")) - 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("Calcite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
    Dolomite
-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("Dolomite") < 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("Dolomite")) - 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("Dolomite")) - 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("Dolomite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Quartz
-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("Quartz") < 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("Quartz")) - 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("Quartz")) - 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("Quartz")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
K-Feldspar
-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("K-Feldspar") < 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("K-Feldspar")) - 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("K-Feldspar")) - 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("K-Feldspar")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Siderite
-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("Siderite") < 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("Siderite")) - 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("Siderite")) - 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("Siderite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Pyrite
-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("Pyrite") < 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("Pyrite")) - 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("Pyrite")) - 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("Pyrite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 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
Illite
-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("Illite") < 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("Illite")) - 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("Illite")) - 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("Illite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Albite
-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("Albite") < 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("Albite")) - 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("Albite")) - 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("Albite")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Smectite-high-Fe-Mg
-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("Smectite-high-Fe-Mg") < 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("Smectite-high-Fe-Mg")) - 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("Smectite-high-Fe-Mg")) - 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("Smectite-high-Fe-Mg")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
Gypsum
-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("Gypsum") < 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("Gypsum")) - 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("Gypsum")) - 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("Gypsum")) - 1)
300 moles = -rate * TIME
310 REM PRINT "kfsp",moles
320 SAVE moles
-end
END
TITLE CO2-H2S Rock Fluid Interaction
SOLUTION 1
    temp      140
    pH        5.5  # Slightly acidic to drive carbonate dissolution
    pe        4
    redox     pe
    units     mg/l
    C(4)      1000  # Higher initial CO2 to prompt dissolution
    Ca        54.6
    Cl        18927
    Fe        0.02
    K         111
    Mg        12.1
    Na        8504
    S(6)      1300
    Sr        0.05
    -water    1 # kg
GAS_PHASE 1
    -fixed_volume               # Switch to a fixed volume to allow CO? to dissolve
    -volume 1                   # Volume in liters
    CO2(g)    81.5292           # CO? partial pressure in atm
    H2S(g)    0.126             # H?S partial pressure in atm
    -temperature 140            # Temperature in ?C

EQUILIBRIUM_PHASES 1
    Albite                   0 190
    Gypsum                   0 100
    Calcite                  0 78200
    Dolomite                 0 60900
    Illite                   0 14400
    K-feldspar               0 152
    Kaolinite                0 84100
    Pyrite                   0 60800
    Quartz                   0 648400
    Siderite                 0 46600
    Smectite-high-Fe-Mg      0 54.6
SAVE SOLUTION 2
END
USE SOLUTION 2

GAS_PHASE 1
    -fixed_pressure
    -pressure 81.6552
    -volume 1
    -temperature 140
    CO2(g)    810.5292
    H2S(g)    0.126

KINETICS 1
Calcite
    -formula  CaCO3  1
    -m        78.2
    -m0       78.2
    -parms    67.1 0.501 1.55e-06 0.0003311 14400 23500 35400 1 1 6 8
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        60.9
    -m0       60.9
    -parms    63.6 6.457e-06 2.9512e-08 7.76e-06 36100 52200 34800 1 1 6 8
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        648.4
    -m0       648.4
    -parms    68.6 0 1.023e-14 0 0 87700 0 0 0 6 8
    -tol      1e-08
K-feldspar
    -formula  KAlSi3O8  1
    -m        0.152
    -m0       0.152
    -parms    71.1 8.71e-11 3.89e-13 6.31e-22 51700 38000 94100 0.5 0.82 6 8
    -tol      1e-08
Siderite
    -formula  FeCO3  1
    -m        46.6
    -m0       46.6
    -parms    46.2 0.000158 1.26e-09 0 45000 62800 0 0.9 0 6 8
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        60.8
    -m0       60.8
    -parms    36.3 4.89e-12 6.6e-14 8.91e-18 65900 22200 17900 0.78 0 6 8
    -tol      1e-08
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
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        14.4
    -m0       14.4
    -parms    460 1.95e-13 3.9e-15 3.9e-15 48000 48000 48000 0.22 -0.13 6 8
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.19
    -m0       0.19
    -parms    69.4 6.91e-11 2.75e-13 0 65000 69800 0 0.46 0 6 8
    -tol      1e-08
Smectite-high-Fe-Mg
    -formula  Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12  1
    -m        0.0546
    -m0       0.0546
    -parms    1040 1.047e-11 1.659e-13 3.019e-17 23600 0 58900 0.34 -0.4 6 8
    -tol      1e-08
Gypsum
    -formula  CaSO4:2H2O  1
    -m        0.1
    -m0       0.1
    -parms    61.6 0 0.001622 0 0 0 0 0 0 6 8
    -tol      1e-08
-steps       1e0 1e1 1e2 1e3 1e4 1e5 1e6 1e7 6.3072e7 # Seconds up to exactly 2 year
#INCREMENTAL_REACTIONS True

USER_GRAPH 3
    -headings               time Volume_Change_Calcite Volume_change_Dolomite Volume_change_Quartz \
                            Volume_change_K-feldspar Volume_change_Siderite Volume_change_Pyrite \
                            Volume_change_Kaolinite Volume_change_Illite Volume_change_Albite \
                            Volume_change_Smectite-high-Fe-Mg Volume_change_Gypsum
    -axis_titles            "Time (Days)" "Phase Volume Changes (cm?/100g)"
    -chart_title            "Phase Volume Changes for Each Mineral Over Time"
    -axis_scale x_axis      auto auto auto auto log
-start
10 GRAPH_X TOTAL_TIME / 86400
20 v_calcite = KIN_DELTA("Calcite") * 36.89
30 v_dolomite = KIN_DELTA("Dolomite") * 64.5
40 v_quartz = KIN_DELTA("Quartz") * 22.69
50 v_kfeldspar = KIN_DELTA("K-feldspar") * 108.15
60 v_siderite = KIN_DELTA("Siderite") * 46.2
70 v_pyrite = KIN_DELTA("Pyrite") * 36.3
80 v_kaolinite = KIN_DELTA("Kaolinite") * 99.35
90 v_illite = KIN_DELTA("Illite") * 102.0
100 v_albite = KIN_DELTA("Albite") * 100.2
110 v_smectite = KIN_DELTA("Smectite-high-Fe-Mg") * 1040
120 v_gypsum = KIN_DELTA("Gypsum") * 58.3
130 GRAPH_Y v_calcite, v_dolomite, v_quartz, v_kfeldspar, v_siderite, v_pyrite, v_kaolinite, v_illite, v_albite, v_smectite, v_gypsum
-end
-active true


[ahmalisha]

Thank you so much for your response. Is there any recommendation that I should make to the code to show the volume change ( like increase the reaction time or the volume of CO2 involved)? much appreciation 

[dlparkhurst]

You can experiment.

[Mira]

Hello,
I trying to simulate the interaction between a Shale and different brine condition at different temperature and pressure. The shale consists minerals like Quartz, calcite, dolomite, illite, kaolinite and pyrite. I want to find out the amount of dissolution and precipitation of minerals over time at above variables. I have written the code but i think something is wrong with it. there are some variation in amount of minerals (not very much) in the end but I am not able to see the variation over time.
I am posting my code can you help me with the possible errors?

Code: [Select]

SOLUTION 1
    temp      45
    pH        6
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Cl        0
    Na        0
    -water    1 # kg

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
    Quartz
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=0
 80 E1=0
 90 n1=0
100 rem neutral solution parameters
110 a2=1.98
120 E2=77000
130 rem base solution parameters
140 a3=1.97E+04
150 E3=80000
160 n2=0.34
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("quartz")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
280 moles= rate*Time/3600
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
-end
    Calcite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be foundin the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid soltion parameters
 70 a1=0
 80 E1=0
 90 n1=0
100 rem neutral solution parameters
110 a2=6.59E+03
120 E2=66000
130 rem co2 solution parameters
140 a3=1.04E+07
150 E3=67000
160 n2=1.6
170 rem rate = 0 if no minerals and undersaturated
180 SR_mineral=SR("calcite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1= a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*ACT("HCO3-")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-SR_mineral)*SA*parm(2)
280 moles= Rate*Time/3600
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
-end
    Dolomite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=3.21E+04
 80 E1=46000
 90 n1=0.500
100 rem neutral solution parameters
110 a2=2.97E-03
120 E2=31000
130 rem CO2 solution parameters
140 a3=0
150 E3=0
160 n2=0
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("dolomite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*ACT("HCO3-")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
280 moles= rate*Time/3600
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
-end
    Illite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals
 50 rem M0 is the initial moles of minerals
 60 rem parm(2) is a scaling factor
 70 rem acid solution parameters
 80 a1=1.00E-02
 90 E1=58000
100 n1=0.55
110 rem neutral solution parameters
120 a2=2.00E-05
130 E2=54000
140 rem base solution parameters
150 a3=1.49E-03
160 E3=77000
170 n2=0.35
180 rem rate=0 if no minerals and undersaturated
190 SR_mineral=SR("Illite")
200 if (M<0) then goto 320
210 if (M=0 and SR_mineral<1) then goto 320
220 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
230 if (SA<=0) then SA=1
240 R=8.31451
250 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
260 Rate2=a2*EXP(-E2/R/TK)
270 Rate3=a3*EXP(-E3/R/TK)*(ACT("OH-"))^n2
280 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
290 moles=Rate*Time/3600
300 rem do not dissolve more minerals than present
310 if (moles>M) then moles=M
320 save moles
-end
    Kaolinite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals
 50 rem M0 is the initial moles of minerals
 60 rem parm(2) is a scaling factor
 70 rem acid solution parameters
 80 a1=2.56E-04
 90 E1=43000
100 n1=0.51
110 rem neutral solution parameters
120 a2=5.00E-08
130 E2=38000
140 rem base solution parameters
150 a3=2.87E-03
160 E3=46000
170 n2=0.58
180 rem rate=0 if no minerals and undersaturated
190 SR_mineral=SR("Kaolinite")
200 if (M<0) then goto 320
210 if (M=0 and SR_mineral<1) then goto 320
220 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
230 if (SA<=0) then SA=1
240 R=8.31451
250 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
260 Rate2=a2*EXP(-E2/R/TK)
270 Rate3=a3*EXP(-E3/R/TK)*(ACT("OH-"))^n2
280 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
290 moles=Rate*Time/3600
300 rem do not dissolve more minerals than present
310 if (moles>M) then moles=M
320 save moles
-end
    Pyrite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals
 50 rem M0 is the initial moles of minerals
 60 rem parm(2) is a scaling factor
 70 rem acid solution parameters
 80 a1=2.82E+02
 90 E1=56900
100 n1=-0.500
110 n3=0.500
120 rem neutral solution parameters
130 a2=2.64E+05
140 E2=56900
150 n2=0.500
160 rem rate=0 if no minerals and undersaturated
170 SR_mineral=SR("Pyrite")
180 if (M<0) then goto 290
190 if (M=0 and SR_mineral<1) then goto 290
200 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
210 if (SA<=0) then SA=1
220 R=8.31451
230 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1*ACT("Fe+3")^n3
240 Rate2=a2*EXP(-E2/R/TK)*ACT("O2")
250 Rate=(Rate1+Rate2)*(1-Sr_mineral)*SA*parm(2)
260 moles=Rate*Time/3600
270 rem do not dissolve more minerals than present
280 if (moles>M) then moles=M
290 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
Quartz
    SiO2 + 2H2O = H4SiO4
    log_k     -3.98
    delta_h   5.99 kcal
    -analytical_expression 0.41 0 -1309 0 0 0
    -Vm       22.67 cm3/mol
Calcite
    CaCO3 = CO3-2 + Ca+2
    log_k     -8.48
    delta_h   -2.297 kcal
    -analytical_expression 17.118 -0.046528 -3496 0 0 0
    -Vm       36.9 cm3/mol
Dolomite
    CaMg(CO3)2 = Ca+2 + Mg+2 + 2 CO3-2
    log_k     -17.09
    delta_h   -9.436 kcal
    -analytical_expression 31.283 -0.0898 -6438 0 0 0
    -Vm       64.5 cm3/mol
Illite
    K0.6Mg0.25Al2.3Si3.5O10(OH)2 + 11.2H2O = 0.6K+ + 0.25Mg+2 + 2.3Al(OH)4- + 3.5H4SiO4 + 1.2H+
    log_k     -40.267
    delta_h   54.684 kcal
    -Vm       141.48 cm3/mol
Kaolinite
    Al2Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 2 Al+3
    log_k     7.435
    delta_h   -35.3 kcal
    -Vm       99.35 cm3/mol
Pyrite
    FeS2 + 2 H+ + 2 e- = Fe+2 + 2 HS-
    log_k     -18.479
    delta_h   11.3 kcal
    -Vm       23.48 cm3/mol

GAS_PHASE 1
    -fixed_volume
    -pressure 70
    -volume 1
    -temperature 45
    CO2(g)    70
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        0.8255
    -m0       0.8255
    -parms    0.0687 1
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.0809
    -m0       0.0809
    -parms    0.0671 1
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0.0141
    -m0       0.0141
    -parms    0.0635 1
    -tol      1e-08
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        0.0217
    -m0       0.0217
    -parms    0.04615 1
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0.019
    -m0       0.019
    -parms    1.1563 1
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        0.01
    -m0       0.01
    -parms    0.0363 1
    -tol      1e-08

-steps       2592000 in 50 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-008
    -tolerance             1e-015
    -step_size             100
    -pe_step_size          10
USER_GRAPH 1
    -axis_titles            "Time" "KIN" ""
    -chart_title            "Dissolution"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TIME/86400
20 GRAPH_Y KIN("Quartz"),KIN("Calcite"),KIN("Dolomite"),KIN("Illite"),KIN("Kaolinite"),KIN("Pyrite")
  -end
    -active                 true

SELECTED_OUTPUT 1
    -file                 C:\Users\hp\Desktop\Mineral_mixture\Shale.sel
    -reset                false
USER_PUNCH 1
    -headings Time(day) pH Mg_conc Ca_conc Quartz Calcite Dolomite Illite Kaolinite Pyrite
    -start
 10 PUNCH KIN("Time")
 20 PUNCH -LA("H+")
 30 PUNCH TOT("Mg")
 40 PUNCH TOT("Ca")
 50 PUNCH KIN("Quartz")
 60 PUNCH KIN("Calcite")
 70 PUNCH KIN("Dolomite")
 80 PUNCH KIN("Illite")
 90 PUNCH KIN("Kaolinite")
100 PUNCH KIN("Pyrite")
    -end
END

 

[dlparkhurst]

Use the Basic variable TOTAL_TIME for your X axis. TIME is an internally generated time step.

Otherwise, you will have to decide if your rates are reacting fast enough for your expectations.

[Mira]

Thank you for the help...now its working.
The dissolution and precipitation of minerals are very minimal...can u suggest me that what could be the reason and how can i get significant change in it over long period of time.

[dlparkhurst]

To get more reaction, consider longer time periods or increase the rates (usually by adjusting surface area).

[Mira]

Hello,
I am trying to simulate the interaction between a basalt and different brine condition for CO2 sequestration.  I have written the code but the code is not running I think something is wrong with it.
I am posting my code can you please help me with the possible errors?

Code: [Select]

RATES
    Albite
-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=1.45E+0
 12 E1=58400
 13 n1=0.335
 20 rem neutral solution parameters
 21 a2=4.97E-10
 22 E2=57000
 30 rem base solution parameters
 31 a3=7.15E-01
 32 E3=55500
 33 n2=0.317
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Albite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
 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
    anorthite
-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.58E-01
 12 E1=16601
 13 n1=1.411
 20 rem neutral solution parameters
 21 a2=1.00E-06
 22 E2=17821
 30 rem base solution parameters
 31 a3=1.00E-22
 32 E3=18150
 33 n2=-1.767
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("anorthite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
    enstatite
-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=1.00E+05
12 E1=80000
13 n1=0.600
20 rem neutral solution parameters
21 a2=2.00E+01
22 E2=80000
30 rem base solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("enstatite")
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("H+")^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
    diopside
-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=3.00E+10
 12 E1=96100
 13 n1=0.710
 20 rem neutral solution parameters
 21 a2=1.00E-04
 22 E2=40600
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("diopside")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
    forsterite
-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=8.38E+04
 12 E1=67206
 13 n1=0.470
 20 rem neutral solution parameters
 21 a2=1.58E+03
 22 E2=79000
 30 rem base solution parameters
 31 a3=1.00E-07
 32 E3=56637
 33 n2=-0.600
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("forsterite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 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
    Magnetite
-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=4.66E-06
 12 E1=18600
 13 n1=0.279
 20 rem neutral solution parameters
 21 a2=3.01E-08
 22 E2=18600
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Magnetite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

  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
    magnesite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=1.39E-4
 80 E1=14400
 90 n1=1.000
100 rem neutral solution parameters
110 a2=5.99E-6
120 E2=23500
130 rem CO2 denpendence parameters
140 a3=6.03E+05
150 E3=62800
160 n2=1.000
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("magnesite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*SR("CO2(g)")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral^4)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
-end

  dolomite
-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=3.21E+04
12 E1=46000
13 n1=0.500
20 rem neutral solution parameters
21 a2=2.97E-03
22 E2=31000
30 rem CO2 solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("dolomite")
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

    Siderite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=3.82E+04
 80 E1=56000
 90 n1=0.6
100 rem neutral solution parameters
110 a2=1.36E+01
120 E2=56000
130 rem base solution parameters
140 a3=0
150 E3=0
160 n2=0
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("Siderite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
-end

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

PHASES
Albite
        NaAlSi3O8 +4.0000 H+  =  + 1.0000 Al+++ + 1.0000 Na+ + 2.0000 H2O + 3.0000 SiO2
        log_k           2.7645
-delta_H -51.8523 kJ/mol # Calculated enthalpy of reaction Albite
# Enthalpy of formation: -939.68 kcal/mol
        -analytic -1.1694e+001 1.4429e-002 1.3784e+004 -7.2866e+000 -1.6136e+006
#       -Range:  0-300

Anorthite
        CaAl2(SiO4)2 +8.0000 H+  =  + 1.0000 Ca++ + 2.0000 Al+++ + 2.0000 SiO2 + 4.0000 H2O
        log_k           26.5780
-delta_H -303.039 kJ/mol # Calculated enthalpy of reaction Anorthite
# Enthalpy of formation: -1007.55 kcal/mol
        -analytic 3.9717e-001 -1.8751e-002 1.4897e+004 -6.3078e+000 -2.3885e+005
#       -Range:  0-300

Enstatite
        MgSiO3 +2.0000 H+  =  + 1.0000 H2O + 1.0000 Mg++ + 1.0000 SiO2
        log_k           11.3269
-delta_H -82.7302 kJ/mol # Calculated enthalpy of reaction Enstatite
# Enthalpy of formation: -369.686 kcal/mol
        -analytic -4.9278e+001 -3.2832e-003 9.5205e+003 1.4437e+001 -5.4324e+005
#       -Range:  0-300

Diopside
        CaMgSi2O6 +4.0000 H+  =  + 1.0000 Ca++ + 1.0000 Mg++ + 2.0000 H2O + 2.0000 SiO2
        log_k           20.9643
-delta_H -133.775 kJ/mol # Calculated enthalpy of reaction Diopside
# Enthalpy of formation: -765.378 kcal/mol
        -analytic 7.1240e+001 1.5514e-002 8.1437e+003 -3.0672e+001 -5.6880e+005
#       -Range:  0-300

Forsterite
        Mg2SiO4 +4.0000 H+  =  + 1.0000 SiO2 + 2.0000 H2O + 2.0000 Mg++
        log_k           27.8626
-delta_H -205.614 kJ/mol # Calculated enthalpy of reaction Forsterite
# Enthalpy of formation: -520 kcal/mol
        -analytic -7.6195e+001 -1.4013e-002 1.4763e+004 2.5090e+001 -3.0379e+005
#       -Range:  0-300

Magnetite
        Fe3O4 +8.0000 H+  =  + 1.0000 Fe++ + 2.0000 Fe+++ + 4.0000 H2O
        log_k           10.4724
-delta_H -216.597 kJ/mol # Calculated enthalpy of reaction Magnetite
# Enthalpy of formation: -267.25 kcal/mol
        -analytic -3.0510e+002 -7.9919e-002 1.8709e+004 1.1178e+002 2.9203e+002
#       -Range:  0-300

Calcite
        CaCO3 +1.0000 H+  =  + 1.0000 Ca++ + 1.0000 HCO3-
        log_k           1.8487
-delta_H -25.7149 kJ/mol # Calculated enthalpy of reaction Calcite
# Enthalpy of formation: -288.552 kcal/mol
        -analytic -1.4978e+002 -4.8370e-002 4.8974e+003 6.0458e+001 7.6464e+001
#       -Range:  0-300

Dolomite
        CaMg(CO3)2 +2.0000 H+  =  + 1.0000 Ca++ + 1.0000 Mg++ + 2.0000 HCO3-
        log_k           2.5135
-delta_H -59.9651 kJ/mol # Calculated enthalpy of reaction Dolomite
# Enthalpy of formation: -556.631 kcal/mol
        -analytic -3.1782e+002 -9.8179e-002 1.0845e+004 1.2657e+002 1.6932e+002
#       -Range:  0-300

Magnesite
        MgCO3 +1.0000 H+  =  + 1.0000 HCO3- + 1.0000 Mg++
        log_k           2.2936
-delta_H -44.4968 kJ/mol # Calculated enthalpy of reaction Magnesite
# Enthalpy of formation: -265.63 kcal/mol
        -analytic -1.6665e+002 -4.9469e-002 6.4344e+003 6.5506e+001 1.0045e+002
#       -Range:  0-300

Siderite
        FeCO3 +1.0000 H+  =  + 1.0000 Fe++ + 1.0000 HCO3-
        log_k           -0.1920
-delta_H -32.5306 kJ/mol # Calculated enthalpy of reaction Siderite
# Enthalpy of formation: -179.173 kcal/mol
        -analytic -1.5990e+002 -4.9361e-002 5.4947e+003 6.3032e+001 8.5787e+001
#       -Range:  0-300

SOLUTION 1
-temp 50
pH   7.5
Na  2.91
Ca  0.521
K   0.039
F  0.017
Cl  2.6
S   0.78  as  SO4
C 4.74 as HCO3-
Si   0.67   as   SiO2

GAS_PHASE 1
    -fixed_pressure
    -pressure 25
    -volume 1
    -temperature 50
    CO(g)     25
KINETICS 1
Albite
    -formula  NaAlSi3O8  1
    -m        0.633
    -m0       0.633
    -parms    4.8 1
    -tol      1e-08
anorthite
    -formula  anorthite  1
    -m        0.647
    -m0       0.647
    -parms    5.3 1
    -tol      1e-08
enstatite
    -formula  enstatite  1
    -m        1.55166
    -m0      1.55166
    -parms    2.3 1
    -tol      1e-08
diopside
    -formula  diopside  1
    -m        0.953
    -m0       0.953
    -parms    1.92 1
    -tol      1e-08
forsterite
    -formula  forsterite  1
    -m        0.457
    -m0       0.457
    -parms    2.1 1
    -tol      1e-08
Magnetite
    -formula  Fe3O4  1
    -m        0.0862
    -m0       0.0862
    -parms    1.5 1
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
Magnesite
    -formula  MgCO3  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
Siderite
    -formula  FeCO3  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
-steps       3000000000 in 20 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
 INCREMENTAL_REACTIONS true
SELECTED_OUTPUT 1
    -file                 C:\Users\hp\Desktop\Mineral_mixture\basalt\Basalt_.sel
    -reset                false
USER_PUNCH 1
    -headings TIME(Years) KIN("Albite") KIN("anorthite") KIN("enstatite") KIN("diopside") KIN("forsterite") KIN("Magnetite") KIN("Calcite") KIN("Dolomite") KIN("Magnesite") KIN("Siderite")
    -start
 10 PUNCH total_TIME/31536000
 20 PUNCH KIN("Albite")*1000
 30 PUNCH KIN("anorthite")*1000
 40 PUNCH KIN("enstatite")*1000
 50 PUNCH KIN("diopside")*1000
 60 PUNCH KIN("forsterite")*1000
 70 PUNCH KIN("Magnetite")*1000
 80 PUNCH KIN("Calcite")*1000
 90 PUNCH KIN("Dolomite")*1000
100 PUNCH KIN("Magnesite")*1000
110 PUNCH KIN("Siderite")*1000
    -end
INCREMENTAL_REACTIONS true
USER_GRAPH 1
    -headings               rxn KIN("Albite") KIN("anorthite") KIN("enstatite") KIN("diopside") KIN("forsterite") KIN("Magnetite") KIN("Calcite") KIN("Dolomite") KIN("Magnesite") KIN("Siderite")
    -axis_titles            "TIME (year)" "milli moles" ""
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/(3600*24*365)
20 GRAPH_Y KIN("Albite")*1000, KIN("anorthite")*1000,KIN("enstatite")*1000, KIN("diopside")*1000, KIN("forsterite")*1000 KIN("Magnetite") KIN("Calcite")*1000, KIN("Dolomite")*1000, KIN("Magnesite")*1000, KIN("Siderite")*1000
  -end
    -active                 true
 END

[Mira]

Anyone please help me with this

[dlparkhurst]

No whining.

You have defined CO(g) in the gas phase. Assuming you want CO2(g), you probably want to define a large volume so that you do not run out of CO2 during the reaction.

Magnetite comes to equilibrium quickly. The simulation will run much faster if you convert Magnetite to an EQUILIBRIUM_PHASE rather than a KINETIC reaction. The small concentrations of Fe cause the integration to struggle if it is a kinetic reactant.

There are some transitions that require smaller time steps. The following uses a progressive set of time increments to get past the difficult spots.
I don't know what database you are using, but this script runs with llnl.dat.

Code: [Select]

RATES
    Albite
-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=1.45E+0
 12 E1=58400
 13 n1=0.335
 20 rem neutral solution parameters
 21 a2=4.97E-10
 22 E2=57000
 30 rem base solution parameters
 31 a3=7.15E-01
 32 E3=55500
 33 n2=0.317
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Albite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
 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
    anorthite
-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.58E-01
 12 E1=16601
 13 n1=1.411
 20 rem neutral solution parameters
 21 a2=1.00E-06
 22 E2=17821
 30 rem base solution parameters
 31 a3=1.00E-22
 32 E3=18150
 33 n2=-1.767
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("anorthite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
    enstatite
-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=1.00E+05
12 E1=80000
13 n1=0.600
20 rem neutral solution parameters
21 a2=2.00E+01
22 E2=80000
30 rem base solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("enstatite")
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("H+")^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
    diopside
-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=3.00E+10
 12 E1=96100
 13 n1=0.710
 20 rem neutral solution parameters
 21 a2=1.00E-04
 22 E2=40600
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("diopside")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
    forsterite
-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=8.38E+04
 12 E1=67206
 13 n1=0.470
 20 rem neutral solution parameters
 21 a2=1.58E+03
 22 E2=79000
 30 rem base solution parameters
 31 a3=1.00E-07
 32 E3=56637
 33 n2=-0.600
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("forsterite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 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
    Magnetite
-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=4.66E-06
 12 E1=18600
 13 n1=0.279
 20 rem neutral solution parameters
 21 a2=3.01E-08
 22 E2=18600
 30 rem base solution parameters
 31 a3=0
 32 E3=0
 33 n2=0
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("Magnetite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

  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
    magnesite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=1.39E-4
 80 E1=14400
 90 n1=1.000
100 rem neutral solution parameters
110 a2=5.99E-6
120 E2=23500
130 rem CO2 denpendence parameters
140 a3=6.03E+05
150 E3=62800
160 n2=1.000
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("magnesite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*SR("CO2(g)")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral^4)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
#300 if (moles>M) then moles=M
310 save moles
-end

  dolomite
-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=3.21E+04
12 E1=46000
13 n1=0.500
20 rem neutral solution parameters
21 a2=2.97E-03
22 E2=31000
30 rem CO2 solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("dolomite")
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

    Siderite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=3.82E+04
 80 E1=56000
 90 n1=0.6
100 rem neutral solution parameters
110 a2=1.36E+01
120 E2=56000
130 rem base solution parameters
140 a3=0
150 E3=0
160 n2=0
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("Siderite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
#300 if (moles>M) then moles=M
310 save moles
-end

SOLUTION_SPECIES
#Al+3 + 4H2O = Al(OH)4- + 4H+
#log_k   -22.7
#H4SiO4 = H4SiO4
# -dw  1.10e-9
# -Vm  10.5  1.7  20  -2.7  0.1291 # supcrt + 2*H2O in a1
SiO2 + 2H2O = H4SiO4
        log_k   -2.71
        delta_h 3.91    kcal
2Al+3 + 2H4SiO4 + H2O = Al2Si2O5(OH)4 + 6H+

PHASES
Albite
        NaAlSi3O8 +4.0000 H+  =  + 1.0000 Al+++ + 1.0000 Na+ + 2.0000 H2O + 3.0000 SiO2
        log_k           2.7645
-delta_H -51.8523 kJ/mol # Calculated enthalpy of reaction Albite
# Enthalpy of formation: -939.68 kcal/mol
        -analytic -1.1694e+001 1.4429e-002 1.3784e+004 -7.2866e+000 -1.6136e+006
#       -Range:  0-300

Anorthite
        CaAl2(SiO4)2 +8.0000 H+  =  + 1.0000 Ca++ + 2.0000 Al+++ + 2.0000 SiO2 + 4.0000 H2O
        log_k           26.5780
-delta_H -303.039 kJ/mol # Calculated enthalpy of reaction Anorthite
# Enthalpy of formation: -1007.55 kcal/mol
        -analytic 3.9717e-001 -1.8751e-002 1.4897e+004 -6.3078e+000 -2.3885e+005
#       -Range:  0-300

Enstatite
        MgSiO3 +2.0000 H+  =  + 1.0000 H2O + 1.0000 Mg++ + 1.0000 SiO2
        log_k           11.3269
-delta_H -82.7302 kJ/mol # Calculated enthalpy of reaction Enstatite
# Enthalpy of formation: -369.686 kcal/mol
        -analytic -4.9278e+001 -3.2832e-003 9.5205e+003 1.4437e+001 -5.4324e+005
#       -Range:  0-300

Diopside
        CaMgSi2O6 +4.0000 H+  =  + 1.0000 Ca++ + 1.0000 Mg++ + 2.0000 H2O + 2.0000 SiO2
        log_k           20.9643
-delta_H -133.775 kJ/mol # Calculated enthalpy of reaction Diopside
# Enthalpy of formation: -765.378 kcal/mol
        -analytic 7.1240e+001 1.5514e-002 8.1437e+003 -3.0672e+001 -5.6880e+005
#       -Range:  0-300

Forsterite
        Mg2SiO4 +4.0000 H+  =  + 1.0000 SiO2 + 2.0000 H2O + 2.0000 Mg++
        log_k           27.8626
-delta_H -205.614 kJ/mol # Calculated enthalpy of reaction Forsterite
# Enthalpy of formation: -520 kcal/mol
        -analytic -7.6195e+001 -1.4013e-002 1.4763e+004 2.5090e+001 -3.0379e+005
#       -Range:  0-300

Magnetite
        Fe3O4 +8.0000 H+  =  + 1.0000 Fe++ + 2.0000 Fe+++ + 4.0000 H2O
        log_k           10.4724
-delta_H -216.597 kJ/mol # Calculated enthalpy of reaction Magnetite
# Enthalpy of formation: -267.25 kcal/mol
        -analytic -3.0510e+002 -7.9919e-002 1.8709e+004 1.1178e+002 2.9203e+002
#       -Range:  0-300

Calcite
        CaCO3 +1.0000 H+  =  + 1.0000 Ca++ + 1.0000 HCO3-
        log_k           1.8487
-delta_H -25.7149 kJ/mol # Calculated enthalpy of reaction Calcite
# Enthalpy of formation: -288.552 kcal/mol
        -analytic -1.4978e+002 -4.8370e-002 4.8974e+003 6.0458e+001 7.6464e+001
#       -Range:  0-300

Dolomite
        CaMg(CO3)2 +2.0000 H+  =  + 1.0000 Ca++ + 1.0000 Mg++ + 2.0000 HCO3-
        log_k           2.5135
-delta_H -59.9651 kJ/mol # Calculated enthalpy of reaction Dolomite
# Enthalpy of formation: -556.631 kcal/mol
        -analytic -3.1782e+002 -9.8179e-002 1.0845e+004 1.2657e+002 1.6932e+002
#       -Range:  0-300

Magnesite
        MgCO3 +1.0000 H+  =  + 1.0000 HCO3- + 1.0000 Mg++
        log_k           2.2936
-delta_H -44.4968 kJ/mol # Calculated enthalpy of reaction Magnesite
# Enthalpy of formation: -265.63 kcal/mol
        -analytic -1.6665e+002 -4.9469e-002 6.4344e+003 6.5506e+001 1.0045e+002
#       -Range:  0-300

Siderite
        FeCO3 +1.0000 H+  =  + 1.0000 Fe++ + 1.0000 HCO3-
        log_k           -0.1920
-delta_H -32.5306 kJ/mol # Calculated enthalpy of reaction Siderite
# Enthalpy of formation: -179.173 kcal/mol
        -analytic -1.5990e+002 -4.9361e-002 5.4947e+003 6.3032e+001 8.5787e+001
#       -Range:  0-300
END
SOLUTION 1
-temp 50
pH   7.5
Na  2.91
Ca  0.521
K   0.039
F  0.017
Cl  2.6
S   0.78  as  SO4
C 4.74 as HCO3-
Si   0.67   as   SiO2

GAS_PHASE 1
    -fixed_pressure
    -pressure 25
    -volume 1e5
    -temperature 50
    CO2(g)     25
EQUILIBRIUM_PHASES
Magnetite 0 0.1
KINETICS 1
Albite
    -formula  NaAlSi3O8  1
    -m        0.633
    -m0       0.0633
    -parms    4.8 1
    -tol      1e-08
anorthite
    -formula  anorthite  1
    -m        0.647
    -m0       0.0647
    -parms    5.3 1
    -tol      1e-08
enstatite
    -formula  enstatite  1
    -m        1.55166
    -m0       0.155166
    -parms    2.3 1
    -tol      1e-08
diopside
    -formula  diopside  1
    -m        0.953
    -m0       0.0953
    -parms    1.92 1
    -tol      1e-08
forsterite
    -formula  forsterite  1
    -m        0.457
    -m0       0.0457
    -parms    2.1 1
    -tol      1e-08
#Magnetite
#    -formula  Fe3O4  1
#    -m        0.0862
#    -m0       0.00862
#    -parms    1.5 1
#    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
Magnesite
    -formula  MgCO3  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
Siderite
    -formula  FeCO3  1
    -m        0
    -m0       0
    -parms    1 1
    -tol      1e-08
-steps       1e2 1e3 1e4 100*1e5 10e6 10*1e7 10*1e8  2*1e9 #3e9 in 20 steps # seconds
-cvode true
INCREMENTAL_REACTIONS true

USER_GRAPH 2
    -headings               TIME(Years) SI("Albite") SI("anorthite") SI("enstatite") SI("diopside") SI("forsterite") SI("Magnetite") SI("Calcite") SI("Dolomite") SI("Magnesite") SI("Siderite")
    -axis_titles            "Time, years" "SI" ""
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 GRAPH_X total_TIME/31536000
 20 GRAPH_Y SI("Albite")
 30 GRAPH_Y SI("anorthite")
 40 GRAPH_Y SI("enstatite")
 50 GRAPH_Y SI("diopside")
 60 GRAPH_Y SI("forsterite")
 70 GRAPH_Y SI("Magnetite")
 80 GRAPH_Y SI("Calcite")
 90 GRAPH_Y SI("Dolomite")
100 GRAPH_Y SI("Magnesite")
110 GRAPH_Y SI("Siderite")
  -end
    -active                 true

USER_GRAPH 1
    -headings               rxn KIN("Albite") KIN("anorthite") KIN("enstatite") KIN("diopside") KIN("forsterite") KIN("Magnetite") KIN("Calcite") KIN("Dolomite") KIN("Magnesite") KIN("Siderite")
    -axis_titles            "TIME (year)" "milli moles" ""
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/(3600*24*365)
20 GRAPH_Y KIN("Albite")*1000, KIN("anorthite")*1000,KIN("enstatite")*1000, KIN("diopside")*1000, KIN("forsterite")*1000 KIN("Magnetite") KIN("Calcite")*1000, KIN("Dolomite")*1000, KIN("Magnesite")*1000, KIN("Siderite")*1000
  -end
    -active                 true
 END


[Mira]

Thank you, now it is working ...
I want the dissolution of minerals and precipitation of carbonates..what change should i make to obtaine that?
can you give me some suggestions for my objective?

[dlparkhurst]

There is a lot of literature on CO2 sequestration.