amount of CO2 injected

[nurlan]

Dear Sir David,
Hope you are doing well,
I have two questions regarding the Co2 injection RTM model.
1. CO2 amount is transported from cell to cell as CO2(g) or as a C element? I sequester the CO2 using equilibrium_phases. When I do the simulation amount of CO2 is rapidly decreases after 10-20 years, and reaches zero after that. But amount of C is never zero. Does that mean Co2 injected stopped or it is a function of C element? because I am worried that I am injecting Co2 amount only once and whole process just monitoring.
2.
I did the dissolution of CO2 in the brine water first, I saved it as Solution 0. brine itself was Solution 1-n. Then, to support the CO2 partial pressure I added Equilibrium_phases with Co2 fugacity parameter. Did I make it correctly, or I should remove the Equil_Phases?
Thus, code looks as below

Code: [Select]

DATABASE c:\phreeqc\database\llnl.dat
SOLUTION 0 #solution injected into the 1st cell 

temp 40
pH 6.55
pe -3.1
redox pe
units mol/kgw
density 1.1596
pressure 88.8

Al 1.15E-08
C 9.88E-01
Ca 3.86E-02
Cl 5.24E+00
Fe 6.84E-04
K 4.02E-02
Mg 6.32E-02
Na 4.90E+00
S 1.17E-06
Si 1.84E-04
Sr 7.82E-04
   
water 1 

#All rate data from Palandri and Kharaka, except for dawsonite...
 Equilibrium_phases 1
 CO2(g) 1.75

SOLUTION 1-50 #solution into which CO2 is dissolved

temp 40
pH 6.55
pe -3.1
redox pe
units mol/kgw
density 1.1596
pressure 88.8

Al            1.153e-08   
C                 9.941e-03   
Ca                3.857e-02   
Cl                5.237e+00   
Fe                6.842e-04   
K                 4.024e-02   
Mg                6.319e-02   
Na                4.903e+00 
S                 1.170e-06   
Si                1.843e-04   
Sr                7.824e-04   
water 1 
   
   
RATES

Quartz
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.0
15 logkH = -13.99
    16 Ea = 87700
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Quartz") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Quartz")))^2)
    24 k_n = 1
25 gamma_n = 5e10
26  lamda = 1
    30  if SR("Quartz") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Quartz"))*time
    40  if SR("Quartz") > 1 then moles = (-reacS*lamda*k_p*(SR("Quartz") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end



Albite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.457
15 logkH = -10.16
    16 Ea = 65000
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Albite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Albite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Albite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Albite"))*time
    40  if SR("Albite") > 1 then moles = (-reacS*lamda*k_p*(SR("Albite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

Kaolinite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.777
15 logkH = -11.31
    16 Ea = 65900
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/10
23 if SR("Kaolinite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Kaolinite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Kaolinite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Kaolinite"))*time
    40  if SR("Kaolinite") > 1 then moles = (-reacS*lamda*k_p*(SR("Kaolinite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

Dolomite-dis #Same as dolomite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.5
15 logkH = -3.19
    16 Ea = 56700
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = 9.45E-15
23 if SR("Dolomite-dis") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Dolomite-dis")))^2)
    24 k_n = 1
25 gamma_n = 4e10
26  lamda = 1
    30  if SR("Dolomite-dis") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Dolomite-dis"))*time
    40  if SR("Dolomite-dis") > 1 then moles = (-reacS*lamda*k_p*(SR("Dolomite-dis") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

Magnesite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 1.0
15 logkH = -6.38
    16 Ea = 14400
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = 7.36E-14
23 if SR("Magnesite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Magnesite")))^2)
    24 k_n = 1
25 gamma_n = 4e10
26  lamda = 1
    30  if SR("Magnesite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Magnesite"))*time
    40  if SR("Magnesite") > 1 then moles = (-reacS*lamda*k_p*(SR("Magnesite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

Siderite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.94
15 logkH = -7.5
    16 Ea = 48000
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Siderite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Siderite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Siderite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Siderite"))*time
    40  if SR("Siderite") > 1 then moles = (-reacS*lamda*k_p*(SR("Siderite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
Calcite #Carbonate mechanism instead of base!!!
-start
10 reacS = parm(1)*parm(2)*m*0.1
11 n_acid = 1
12 n_carb = 1
13 logkH = -0.30
14 logkH2O = -5.81
15 logkCARB = -3.48
        16 EaH = 14400
17 EaH2O = 23500
18 EaCARB = 35400
        19 R = 8.3145
        20 logkHT = logkH - (EaH/(2.303*R))*(1/Tk - 1/298.15)
        21 logkH2OT = logkH2O - (EaH2O/(2.303*R))*(1/Tk - 1/298.15)
22 logkCARBT = logkCARB - (EaCARB/(2.303*R))*(1/Tk - 1/298.15) #
23 kHT = 10^logkHT
        24 kH2OT = 10^logkH2OT
25 kCARBT = 10^logkCARBT
26 k_p = ((act("H+")^n_acid)*kHT + kH2OT + (act("HCO3-")^n_carb)*kCARBT)
        30  if SR("Calcite") <= 1 then moles = reacS*k_p*(1 - SR("Calcite"))*time
        40  if SR("Calcite") > 1 then moles = (-(reacS + 0.001)*k_p/100*(SR("Calcite") - 1)^2)*time
        100 save moles
-end
Dawsonite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.98
15 logkH = -4.5
    16 Ea = 63820
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT
23 if SR("Dawsonite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Dawsonite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Dawsonite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Dawsonite"))*time
    40  if SR("Dawsonite") > 1 then moles = (-reacS*lamda*k_p*(SR("Dawsonite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

KINETICS 1-50
-cvode
Quartz
-m0    158.2
-parms  60.08  58


Albite
-m0    9.0
-parms 263.02  21.6



Kaolinite
-m0    3.1
-parms  258.16   200


Dolomite-dis
  -m0     0.0
-parms  206.336    2.8

Magnesite
          -m0     0.177
          -parms   84.31  3.85
Calcite           
-m0 21.90
-parms  100.00  22

Siderite
-m0     0.0
  -parms 115.86     0.175



Dawsonite
-m0     0.0
        -parms 144.0  9.8

            -bad_step_max     2000
           
PHASES

#Valid for 37C
Dawsonite
        NaAlCO3(OH)2 +3.0000 H+  =  + 1.0000 Al+++ + 1.0000 HCO3- + 1.0000 Na+ + 2.0000 H2O
        log_k           4.3464
-delta_H -76.3549 kJ/mol

        -analytic -1.1393e+002 -2.3487e-002 7.1758e+003 4.0900e+001 1.2189e+002

TRANSPORT
-cells 50 # total length = cells X length
-length 1.0 # total depth = 100x1.0=100m
-shifts 400 # =
-time_step 3.1536e+8# 20 years in seconds
-flow_direction forward
-boundary_conditions flux flux
-diffusion_coefficient  3e-9
-dispersivities 0.1 # (m)
-correct_disp true
-punch_cells 1 25 50
-punch_frequency 1
-print_cells 1-50
-print_frequency 1

SELECTED_OUTPUT
-file Kz_RTM_Eb-1 So Eq diff.sel
-pe true
-totals Ca Mg Na K Si C(4) Al Fe O(0) Fe(3) Fe(2) S(6)
-si CO2(g) Quartz Albite Kaolinite  Magnesite Dawsonite Dolomite-dis
-kinetic_reactants  Quartz Albite Calcite Kaolinite Magnesite Dawsonite Siderite Dolomite-dis
-equilibrium_phases CO2(g)

INCREMENTAL_REACTIONS true

END
Thanks for your help and wish a good day
 

[dlparkhurst]

SOLUTION 0 has about 1 mol/kgw of carbon (and ~1 kg water), so 1 mol of carbon (mostly CO2) is advected into the column per shift throughout the simulation.

EQUILIBRIUM_PHASES 1 has 10 moles of CO2 to react. It will dissolve (or exsolve) as necessary to maintain a partial pressure of 10^-1.75 in cell 1. If 10 moles dissolve, then there is no more CO2 available to react in cell 1.

No gas transport occurs. All carbon transport is through the dissolved species. There will always be some dissolved carbon in solution, depending on the rates of reaction and ultimately, equilibrium with the minerals that are being reacted kinetically.

[nurlan]

Thank you, Mr. Parkhurst.
I got confused at first, now it's more clear.