CO2-Calcite-Water interaction -CO2 storage in an aquifer

[MMMM]

Hello
I'm attempting to model a CO2 storage system in an aquifer with 44.2 mol calcite, 1 kg water, and 0.945 mol CO2 at 33 ?C and 987 atm pressure. I'd like to know how much calcite and CO2 are dissolved in water after 240 years.
I wrote the following script but am not sure if it is correct.

I will appreciate anyone who can assist me.

Code: [Select]

[DATABASE C:\phreeqc\database\LLNL.dat
SOLUTION  1      Freshwater
end
EQUILIBRIUM_PHASES  1
Calcite     0 44.2     
Save Solution 1
Save Equilibrium_Phases 1

Use Solution 1
Use Equilibrium_Phases 1

GAS_PHASE 1
     -fixed_pressure
     -pressure 1
     -temperature 33.33
CO2(g) 0
CH4(g) 0
H2O(g) 0
REACTION 1
 CO2(g) 1
0.94522
USE REACTION_PRESSURE 1
USE GAS_PHASE 1
Use Solution 1
Use Equilibrium_Phases 1
REACTION_TEMPERATURE 1
33.33
REACTION_PRESSURE 1
987
#########################################
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


CO2(g) #_hydroxylation
-start
# CO2 + H2O = HCO3- + H+
10 kf = 8.5e-3 # mol-1 s-1
20 kb = 2e-4  # 1/s
30 rf = kf*ACT("[CO2]")*ACT("H2O")
40 rb = kb*ACT("HCO3-")*ACT("H+")
50 rate = (rf - rb)
60 moles = rate * TIME
70 SAVE moles
-end
 

########################################
KINETICS 1
Calcite
     -M0  44.2                   # initial moles   
     -tol 1e-8
     -parms 5 0.3
     

     
CO2(g)  #CO2_hydroxylation
     -M0  0.94522                  # initial moles   
     -tol 1e-8
     -formula CO2 -1 CO2 +1
     
     
 -time  240 year in 24
#############################

USER_GRAPH 5
-headings Title "Dissolved-Calcite"
-axis_titles "Time (year)" "Amount of dissolved minerals, percent of initial"
-axis_scale x_axis auto auto auto aut
-axis_scale y_axis auto auto auto aut
initial_solutions true
-start
10 x = TOTAL_TIME/3.1536e+7
20 GRAPH_X x
30 GRAPH_Y ((EQUI_DELTA("Calcite")*(-1))/44.2)*100
-end

USER_GRAPH 6
-headings CO2
-axis_titles "Time (year)" "CO2 free gas Loss, percent of inventory"
-axis_scale x_axis auto auto auto aut
-axis_scale y_axis auto auto auto aut
initial_solutions true
-start
10 x = TOTAL_TIME/3.1536e+7
20 GRAPH_X x
30 GRAPH_Y ((0.94522-GAS("CO2(g)"))/0.94522)*100
-end

USER_GRAPH 7
-headings Title "Dissolved-CO3--" "Amount-of-dissolved-H+"
-axis_titles "Time (year)" "Amount of dissolved CO3--, gmole" "Amount of dissolved H+, gmole"
-axis_scale x_axis auto auto auto aut
-axis_scale y_axis auto auto auto aut
initial_solutions true
-start
10 x = TOTAL_TIME/3.1536e+7
20 GRAPH_X x
30 GRAPH_Y tot("C(+4)") tot("H")
-end


USER_GRAPH 8
-headings Title "Dissolved-CO2"
-axis_titles "Time (year)" "Amount of dissolved CO2, gmole"
-axis_scale x_axis auto auto auto aut
-axis_scale y_axis auto auto auto aut
initial_solutions true
-start
10 x = TOTAL_TIME/3.1536e+7
20 GRAPH_X x
30 GRAPH_Y mol("CO2")
-end
End
/code]

[dlparkhurst]

Calcite reacts quickly, hours to days at 25 C and 1 atm. I think you can safely assume equilibrium after 240 years, so remove the KINETICS definitions. Also, at least compare results from llnl.dat with phreeqc.dat.