A question about the kinetics simulation of CO2 transport after CO2 injection

[amueru]

I want to do a simulation to see how will the Calcite(the mineral of carbonate reservoir) and CO2 change with time after the injection of CO2 mixed gas.
I write the code like this:

Code: [Select]

SOLUTION 1
ph       7 charge
-units ppm
-pressure    296.07698
-temperature    100
Na       23867
Ca       13961
Mg       2177
Cl       57606
S(6)       5571 as SO4
K       0
C       243 as HCO3
-water        150

GAS_PHASE 1
-fixed_pressure
-pressure 298
-volume   150 #Gas total volume
-temp     100
CO2(g)  127.5 #Litres   
H2S(g)  7.5
SO2(g)  7.5
NO2(g)  7.5
SAVE SOLUTION  1
SAVE GAS_PHASE 1
END

USE SOLUTION 1
KINETICS 1
Calcite
     -M0  148.9805359                  # initial moles   
     -tol 1e-8         
     
Dolomite
     -M0  0                   
     -tol 1e-8         
   
Quartz
     -M0  0               
     -tol 1e-8         
     
Siderite
     -M0  0                   
     -tol 1e-8

#Pyrite
     #-M0  0                   
    #-tol 1e-8

#Anhydrite
    #-M0  0                   
    #-tol 1e-8     
     
-time_step 36500 day in 100             # 100 years
#-step_divide 1000
-cvode  true

EQUILIBRIUM_PHASES 1
#Pyrite 0 0
#Anhydrite 0 0
RATES
Calcite
-start
10 mole = 0
20 If (m <= 0) and (SR("Calcite") < 1) Then GoTo 250
30 S = 0.7 # average BET; suggested value in m2/g
40 Mm = 100.087 # molar mass in g/mol
50 If (SR("Calcite") > 1) Then GoTo 130
########## start dissolution bloc ##########
60 knu = 1.6E-6 * exp((-24000 / 8.314) * ((1 / TK) - (1 / 298.15)))
70 k1 = 5E-1 * exp((-14000 / 8.314) * ((1 / TK) - (1 / 298.15))) * (ACT("H+") ^ 1)
80 k = knu + k1
# kinetic data extracted from 04pal/kha
90 theta = 1 # default value
100 eta = 1 # default value
110 rate = S * m * Mm * k * ((1 - SR("Calcite") ^ theta) ^ eta)
120 GoTo 240
########## end dissolution bloc ##########
########## start precipitation bloc ##########
130 knu = 1.8E-7 * exp((-66000 / 8.314) * ((1 / TK) - (1 / 298.15)))
140 k1 = 1.9E-3 * exp((-67000 / 8.314) * ((1 / TK) - (1 / 298.15))) * (ACT("HCO3-") ^ 1.63)
150 kpre = (-1) * (knu + k1)
160 theta = 0.50
170 eta =  2.00
# kinetic data extracted from 07lio/par 98zha/gra 85ins/blo 95shi/bra 92rom/gro 98zah/daw 98leb/sua
180 If (m <= 0) then GoTo 210
190 rate = S * m * Mm * kpre * (ABS(1 - SR("Calcite") ^ theta) ^ eta)
200 GoTo 240
#start nucleation
210 rate = -1E-10
#end nucleation
########## end precipitation bloc ##########
240 mole = rate * Time
250 Save mole
-end

Dolomite
-start
10 mole = 0
20 If (m <= 0) and (SR("Dolomite") < 1) Then GoTo 240
30 S = 0.09 # average BET; suggested value in m2/g
40 Mm = 184.401 # molar mass in g/mol
50 If (SR("Dolomite") > 1) Then GoTo 130
########## start dissolution bloc ##########
60 knu = 1.1E-8 * exp((-31000 / 8.314) * ((1 / TK) - (1 / 298.15)))
70 k1 = 2.8E-4 * exp((-46000 / 8.314) * ((1 / TK) - (1 / 298.15))) * (ACT("H+") ^ 0.61)
80 k = knu + k1
# kinetic data extracted from 05pok/gol 07gau/sch 99gau/oel 09pok/gol 82bus/plu
90 theta = 0.16 # default value
100 eta = 2.1 # default value
# theta = 0.16 and eta = 2.1 at 80ｰC and pH 7 (extracted from 07gau/sch)
110 rate = S * m * Mm * k * ((1 - SR("Dolomite") ^ theta) ^ eta)
120 GoTo 230
########## end dissolution bloc ##########
########## start precipitation bloc ##########
130 knu = 9.5E-15 * exp((-103000 / 8.314) * ((1 / TK) - (1 / 298.15)))
140 kpre = (-1) * knu
150 theta = 1
160 eta =  1
# kinetic data extracted from 97arv/mac 99arv/mac
170 If (m <= 0) then GoTo 200
180 rate = S * m * Mm * kpre * (ABS(1 - SR("Dolomite") ^ theta) ^ eta)
190 GoTo 230
#start nucleation
200 rate = -1E-10
#end nucleation
########## end precipitation bloc ##########
230 mole = rate * Time
240 Save mole
-end

Quartz
-start
10 mole = 0
20 If (m <= 0) and (SR("Quartz") < 1) Then GoTo 240
30 S = 0.05 # average BET; suggested value in m2/g
40 Mm = 60.084 # molar mass in g/mol
50 If (SR("Quartz") > 1) Then GoTo 130
########## start dissolution bloc ##########
60 knu = 6.4E-14 * exp((-77000 / 8.314) * ((1 / TK) - (1 / 298.15)))
70 k1 = 1.9E-10 * exp((-80000 / 8.314) * ((1 / TK) - (1 / 298.15))) * (ACT("OH-") ^ 0.34)
80 k = knu + k1
# kinetic data extracted from 95kna/cop 88kna/wol 06bic/nag 91ben 88ben/mel 90cas/las 87sch/wal 90bra/wal 00ice/dov 90dov/cre 99dov 94dov 92hou/orr 90blu/yun
90 theta = 1 # default value
100 eta = 1 # default value
110 rate = S * m * Mm * k * ((1 - SR("Quartz") ^ theta) ^ eta)
120 GoTo 230
########## end dissolution bloc ##########
########## start precipitation bloc ##########
130 knu = 3.2E-12 * exp((-50000 / 8.314) * ((1 / TK) - (1 / 298.15)))
140 kpre = (-1) * knu
150 theta = 4.58
160 eta =  0.54
# kinetic data extracted from 05gan/hus 80rim/bar
170 If (m <= 0) then GoTo 200
180 rate = S * m * Mm * kpre * (ABS(1 - SR("Quartz") ^ theta) ^ eta)
190 GoTo 230
#start nucleation
200 rate = -1E-10
#end nucleation
########## end precipitation bloc ##########
230 mole = rate * Time
240 Save mole
-end

Siderite
-start
10 mole = 0
20 If (m <= 0) and (SR("Siderite") < 1) Then GoTo 240
30 S = 2.7 # average BET; suggested value in m2/g
40 Mm = 115.856 # molar mass in g/mol
50 If (SR("Siderite") > 1) Then GoTo 130
########## start dissolution bloc ##########
60 knu = 2.1E-9 * exp((-56000 / 8.314) * ((1 / TK) - (1 / 298.15)))
70 k1 = 5.9E-6 * exp((-56000 / 8.314) * ((1 / TK) - (1 / 298.15))) * (ACT("H+") ^ 0.60)
80 k = knu + k1
# kinetic data extracted from 04duc/mar 11tan/mar 09gol/ben 02pok/sch 86bie/hou 03kop/hom 09tes/duf
90 theta = 1 # default value
100 eta = 1 # default value
110 rate = S * m * Mm * k * ((1 - SR("Siderite") ^ theta) ^ eta)
120 GoTo 230
########## end dissolution bloc ##########
########## start precipitation bloc ##########
130 knu = 1.6E-11 * exp((-108000 / 8.314) * ((1 / TK) - (1 / 298.15)))
140 kpre = (-1) * knu
150 theta = 1
160 eta =  1
# kinetic data extracted from 04jim/rom 92gre/tom
170 If (m <= 0) then GoTo 200
180 rate = S * m * Mm * kpre * (ABS(1 - SR("Siderite") ^ theta) ^ eta)
190 GoTo 230
#start nucleation
200 rate = -1E-10
#end nucleation
########## end precipitation bloc ##########
230 mole = rate * Time
240 Save mole
-end



INCREMENTAL_REACTIONS true
SAVE SOLUTION 2

SELECTED_OUTPUT
        -file          kinetics_results_II.xls
        -temp          true
        -totals        Ca Mg #mol/L
   -molalities    HCO3-       #mol/L
   -kinetics       Calcite Dolomite Quartz Siderite
   -equilibrium_phases  CO2(g) H2S(g) SO2(g) NO2(g) Pyrite Anhydrite
   
PRINT
-selected_out true
END

And I got a result which shows the Calcite will not change with time (k_Calcite shows the same data in each step).
Could you kindly tell me which part(s) I have written is wrong?

[dlparkhurst]

Your calcite rate is sufficiently fast to reach equilibrium in about 10 seconds.

I would not use the additional gases H2S, SO2, and NO2. They will generate some redox reactions that will complicate the interpretation of the results.

[amueru]

Your calcite rate is sufficiently fast to reach equilibrium in about 10 seconds.

I would not use the additional gases H2S, SO2, and NO2. They will generate some redox reactions that will complicate the interpretation of the results.

Thank you for your reply.
I'm trying to do a research about how will the additional gases impact the injection and storage of CO2, so I can't simply delete them.
About the rate, I got them from some former researches. Should I define the rate by myself?

[dlparkhurst]

If you dissolve SO2 or NO2 at equilibrium, they will immediately react to other redox states of S and N. They will not be stable in a solution.

To do any simulations where SO2 concentrations are significant, you will have to separate SO2 and NO2 as separate "elements". That way they are disconnected from other S and N redox states and will be stable. If you need them to react to other redox states, they will have to react through KINETICS or REACTION.

Here is a script that uses data from llnl.dat.  is defined as a new element that has aqueous species of the SO2 redox state. Compare the S and  system in the following calculation.

Code: [Select]

SOLUTION_MASTER_SPECIES
    [S]           [S]O3-2          1     S               12
SOLUTION_SPECIES
[S]O3-2 = [S]O3-2
log_k 0.0

2.0000 H+ + 1.0000 [S]O3-2  =  [S]O2 +1.0000 H2O
        -llnl_gamma           3.0   
        log_k           +9.0656
-delta_H 26.7316 kJ/mol # Calculated enthalpy of reaction SO2
# Enthalpy of formation: -77.194 kcal/mol
        -analytic 9.4048e+001 6.2127e-002 -1.1072e+003 -4.0310e+001 -1.7305e+001
#       -Range:  0-300

2.0000 H+ + 1.0000 [S]O3-2  =  H2[S]O3
        -llnl_gamma           3.0   
        log_k           +9.2132

1.0000 [S]O3-2 + 1.0000 H+  =  H[S]O3-
        -llnl_gamma           4.0   
        log_k           +7.2054
-delta_H 9.33032 kJ/mol # Calculated enthalpy of reaction HSO3-
# Enthalpy of formation: -149.67 kcal/mol
        -analytic 5.5899e+001 3.3623e-002 -5.0120e+002 -2.3040e+001 -7.8373e+000
#       -Range:  0-300

PHASES
                                                                 
[S]O2(g)
[S]O2 = [S]O2
log_k 0.1700
-delta_H 0       # Not possible to calculate enthalpy of reaction SO2(g)
# Enthalpy of formation: 0 kcal/mol
-analytic   -2.0205e+001    2.8861e-003    1.4862e+003    5.2958e+000    1.2721e+005
#       -Range:  0-300
END
SOLUTION
EQUILIBRIUM_PHASES
[S]O2(g) -1 1
SO2(g) -1 1
END


[amueru]

You mean that I should define SO2(g) and "S“Ox for other redox states of S?
Can I use the same data with that I use to define SO2(g)?

And I also have a question about the RATE
I told my problem to the professor and he suggested me to consider about the surface area of each minerals. How can I put the data of surface area into RATE?

[dlparkhurst]

You need to define SO2 as I have done to separate the SO2 species from the rest of the sulfur system. You will have to decide whether you can rely on equilibrium for the remainder of the sulfur system, or if you want to be able to allow more sulfur redox states to exist in non-equilibrium conditions. You will have to define a new "element" for each redox state that you want to exist out of equilibrium. You will then need to use KINETICS or REACTION to allow the individual redox states to react.

For starters, I would probably treat calcite in EQUILIBRIUM_PHASES. If you want a kinetic reaction that considers surface area, look at some of the RATES definitions in phreeqc.dat.

[amueru]

thank you for your reply.
I have already made a research about calcite by equilibrium-phase, and now I want to see the reaction by kinetics.
I tried to use the database of RATE in phreeqc.dat, but I get an error massage

Code: [Select]

ERROR: Parameter subscript out of range.
ERROR:  in BASIC line
 80   IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
ERROR: Fatal Basic error in rate Calcite.
Stopping.

I set the M0 of calcite=148.9805359

[dlparkhurst]

Look at the RATES definition and read the manual.