calcite kinetic rate

[maria]

Dear David,

You had helped me some time ago to set up a synthetic fluid solution for a lab experiment. I am trying now to re-run the file because I have recalculated the mineral moles present in the rock sample that are in contact with the fluid composition and produce a graph in terms of when it will reach equilibrium. Unfortunately, it is not working. I don't understand what I am missing.

Here is my code

Code: [Select]

SOLUTION 1
-pressure 136
-temp 200
SELECTED_OUTPUT 100
-file mg2mol.pqi
USER_PUNCH 100
10 s$ = "REACTION 1" + EOL$
20 s$ = s$ + "  Na2SO4 " + STR_F$(2000*1e-3 / GFW("Na2SO4"), 15, 10) + EOL$
30 s$ = s$ + "  KCl    " + STR_F$(500*1e-3 / GFW("KCl"), 15, 10) + EOL$
40 s$ = s$ + "  CaCO3  " + STR_F$(5.7*1e-3 / GFW("CaCO3"), 15, 10) + EOL$
50 s$ = s$ + "1 mole" + EOL$
60 PUNCH s$
END
SELECTED_OUTPUT 100
-active false
END
USE solution 1
INCLUDE$ mg2mol.pqi
SAVE solution 2
END
MIX
2 0.310
SAVE solution 3
END
USE solution 3
EQUILIBRIUM_PHASES
Calcite 0 24.3
CO2(g) 2.14 1 # adding 1 mol
SAVE solution 4
END

RATES
Calcite
   -start
1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
2   REM   PARM(2) = exponent for M/M0

10  si_cc = SI("Calcite")
20  IF (M <= 0  and si_cc < 0) THEN GOTO 200
30  k1 = 10^(0.198 - 444.0 / TK )
40  k2 = 10^(2.84 - 2177.0 /TK )
50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)
60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )
80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME # convert from mmol to mol
200 SAVE moles
end

INCREMENTAL_REACTIONS
USE solution 4

KINETICS
Calcite
 -m     24.3
 -parms 1e4   0.8  # cm^2/mol calcite, exp factor
 -tol   1e-08
 -time  7200 in 10 steps
USER_GRAPH 1
    -headings               Time SI(Calcite)
    -axis_titles            "Hours" "SI(Calcite)" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y SI("Calcite")
end
Thank you very much in advance.

Kind regards,
Maria

[dlparkhurst]

In your reaction sequence, you have made solution 4 in equilibrium with calcite (EQUILIBRIUM_PHASES). What is the point of a KINETIC reaction? The rate expression is zero at equilibrium.

If you use -cvode in KINETICS, it will run, and show SI ~ 0 for all times. I think the rate is so fast that the Runge-Kutta method of integration has problems; if you run for 1e-3 seconds, it runs successfully with RK, but again shows equilibrium for all times.

[maria]

Dear David,

Thank you very much for your reply. Yes, you are right. I see the point.

But even if I remove calcite under equilibrium phases, and I leave m 24 (as this is the mole amount of calcite present in the rock sample), I am getting some crazy numbers in the y-axis of the SI; SI 0 till 200 x 10-12.

Can this be right?

[dlparkhurst]

Yes, it is right. Your rate is extremely fast, probably because of the large surface area that you have defined. The rate causes calcite equilibrium in less than a fraction of a second, so when you run for hours, the SI is always very close to zero.

Code: [Select]

RATES
Calcite
   -start
1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
2   REM   PARM(2) = exponent for M/M0

10  si_cc = SI("Calcite")
20  IF (M <= 0  and si_cc < 0) THEN GOTO 200
30  k1 = 10^(0.198 - 444.0 / TK )
40  k2 = 10^(2.84 - 2177.0 /TK )
50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)
60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )
80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME # convert from mmol to mol
200 SAVE moles
END
SOLUTION 1
-pressure 136
-temp 200
SELECTED_OUTPUT 100
-file mg2mol.pqi
USER_PUNCH 100
10 s$ = "REACTION 1" + EOL$
20 s$ = s$ + "  Na2SO4 " + STR_F$(2000*1e-3 / GFW("Na2SO4"), 15, 10) + EOL$
30 s$ = s$ + "  KCl    " + STR_F$(500*1e-3 / GFW("KCl"), 15, 10) + EOL$
40 s$ = s$ + "  CaCO3  " + STR_F$(5.7*1e-3 / GFW("CaCO3"), 15, 10) + EOL$
50 s$ = s$ + "1 mole" + EOL$
60 PUNCH s$
END
SELECTED_OUTPUT 100
-active false
END
USE solution 1
INCLUDE$ mg2mol.pqi
SAVE solution 2
END
MIX
2 0.310
SAVE solution 3
END
USE solution 3
EQUILIBRIUM_PHASES
#Calcite 0 24.3
CO2(g) 2.14 1 # adding 1 mol
SAVE solution 4
END



INCREMENTAL_REACTIONS
USE solution 4

KINETICS
Calcite
 -m     24.3
 -parms 1e4   0.8  # cm^2/mol calcite, exp factor
 -tol   1e-08
 #-time  7200 in 10 steps
 -time 3*1e-7 3*1e-6 3*1e-5 3*1e-4 3*1e-3 3*1e-2
 -cvode
USER_GRAPH 1
    -headings               Time SI(Calcite)
    -axis_titles            "Seconds" "SI(Calcite)" ""
    -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")
  -end
    -active                 true
end