long-term behaviour of hydrogen in reservoirs.

[amaury22]

Dear Parkhurst

My goal is to investigate the long-term behavior of hydrogen in reservoirs. In that regard, I am interested in studying the reaction of hydrogen with calcite and H2S, as you suggested we prepare a code. However, I am uncertain about it. By utilizing the forum recommendations and notes, we developed a code that seems to be functioning well and producing satisfactory results. Nevertheless, I am unsure if it is entirely accurate. Could you please review the methodology and provide feedback? Additionally, I would appreciate it if you could guide me to relevant information regarding the reaction between hydrogen Hdg and H2

Code: [Select]

#Kinetic Modeling  Hdg
DATA llnl.dat
#DATABASE llnl.dat

RATES

 Quartz
1 REM d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
2 REM Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
3 REM  k = 10^-13.7 mol/m2/s (25 C)
4 REM  A0, initial surface of quartz (m2) recalc's to mol/s
5 REM  V, solution volume in contact with A0
-start
10  A0 = parm(1)
20  V =  parm(2)
30  dif_temp = 1/TK - 1/298
40  pk_w = 13.7 + 4700.4 * dif_temp
50 rate = (A0 / V) * (m/m0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
60 save rate * time
-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
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
90 theta = 0.16 # default value
100 eta = 2.1 # default value
110 rate = S * m * Mm * k * ((1 - SR("Dolomite") ^ theta) ^ eta)
120 GoTo 230
130 knu = 9.5E-15 * exp((-103000 / 8.314) * ((1 / TK) - (1 / 298.15)))
140 kpre = (-1) * knu
150 theta = 1
160 eta =  1
170 If (m <= 0) then GoTo 200
180 rate = S * m * Mm * kpre * (ABS(1 - SR("Dolomite") ^ theta) ^ eta)
190 GoTo 230
200 rate = -1E-10
230 mole = rate * Time
240 Save mole
-end

K-Feldspar
-start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kfeld = SI("K-Feldspar")
30 if (M <= 0 and si_kfeld < 0) then goto 200             
40 SA = PARM(1) * M             
50 if (M = 0 and si_kfeld > 0) then SA = 1e-05  #nucleation
60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_kfeld))
190 moles = r * TIME
200 SAVE moles
-end

Illite
-start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_ill = SI("Illite")
30 if (M <= 0 and si_ill < 0) then goto 200           
40 SA = PARM(1) * M             
50 if (M = 0 and si_ill > 0) then SA = 1e-05  #nucleation
60 k_acid = 10^(-10.12)*EXP(-58.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("H+")^(0.550)                 
70 k_neut = 10^(-12.26)*EXP(-54.00e+03/8.314*(1.0/TK-1.0/373.15))
80 k_base = 10^(-10.60)*EXP(-77.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("OH-")^(-0.350)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_ill))
190 moles = r * TIME
200 SAVE moles
-end

Kaolinite
-start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kao = SI("Kaolinite")
30 if (M <= 0 and si_kao < 0) then goto 200             
40 SA = PARM(1)*M             
50 if (M = 0 and si_kao > 0) then SA = 1e-05  #nucleation
60 k_acid = 0
70 k_neut = (4.9e-12)*EXP(-65.9e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(si_kao))*ACT("H+")^0.2
190 moles = r * TIME
200 SAVE moles
-end

Chlorite(14A)
-start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_chlo = SI("Chlorite(14A)")
30 if (M <= 0 and si_chlo < 0) then goto 200             
40 SA = PARM(1)*M             
50 if (M = 0 and si_chlo > 0) then SA = 1e-05  #nucleation
60 k_acid = 0
70 k_neut = (7.76e-12)*EXP(-88e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(si_chlo))*ACT("H+")^0.5
190 moles = r * TIME
200 SAVE moles
-end

Albite
-start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alb = SI("Albite")
30 if (M <= 0 and si_alb < 0) then goto 200             
40 SA = PARM(1) * M             
50 if (M = 0 and si_alb > 0) then SA = 1e-05  #nucleation
60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))                 
70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_alb))
190 moles = r * TIME
200 SAVE moles
-end

Calcite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_calc = SI("Calcite")
30 if (M <= 0 and si_calc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_calc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_calc))
190 moles = r * TIME
200 SAVE moles
   -end   


Hdg(g)
-start
10 if (M < 0) then goto 70
20 rate = 2.3e-09*(TOT("C(4)")/1e-3 + TOT("C(4)")) + 9.26e-8*(TOT("S(6)")/(1.e-4 + TOT("S(6)")))
30 moles = rate * TIME
40 if (moles > M) then moles = M
70 SAVE moles
-end


KINETICS

Quartz 1
-formula SiO2
-m0  138                       
-parms 23.13  0.16
-tol 1e-12
 
Dolomite 1
-formula  CaMg(CO3)2
-m0       3.97
-parms    33.192
-tol 1e-12

K-Feldspar 1
-formula  KAlSi3O8  1
-m0       0.07
-parms    86.282
-tol      1e-12

Illite 1
-formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
-m0       0.39
-parms    5956.902
-tol      1e-12

Kaolinite 1
-formula Al2Si2O5(OH)4
-m0 0.072
-parms 298.377
-tol 1e-12

Chlorite(14A) 1
-formula  Mg5Al2Si3O10(OH)8  1
-m0       0.017
-parms    1
-tol      1e-12

Albite 1
-formula  NaAlSi3O8  1
-m0       0.035
-parms    5.260
-tol      1e-12

Calcite
-formula CaCO3 1
    -m 100.09
    -m0 100.09
    -parms 0.001499
    -tol 1e-08

Hdg(g)
     -M0  3.340580489                  # initial moles   
     -tol 1e-8
     -formula Hdg -1 H2 +1

#-time_step 2 day in 2                  # 2 time steps
-cvode   true


-steps 15 year in 300 steps


INCREMENTAL_REACTIONS true
RUN_CELLS
-cell 1


EQUILIBRIUM_PHASES
#Fix_pe 0.78 O2(g)
#CO2(g) 1.9886 #97.4 atm
#Hdg(g)  3
Calcite 0.0 0.185 dissolve only     


SOLUTION 1 Interval Sampled 3-8-22
units mg/L
temp 85.1 #131.1 F
pe -0.78
pH 7.2
Al 0.4
Ba 0.43
Br 9.0
Sr 13.5
Alkalinity 248 as HCO3
Cl 10700
S(6) 2680 charge
F 1.0
Ca 606
Fe(2) 78.5
Mg 110
K 130
Si 11.2
Na 6410

PHASES
CO2(g)
CO2 = CO2
     log_k -1.468
     delta_h -4.776 kcal
     -analytical_expression 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5
     -T_c 304.2 # critical T, K
     -P_c 72.86 # critical P, atm
     -Omega 0.225 # acentric factor
 

EQUILIBRIUM_PHASES 1
CO2(g) 2.10 0.08345 # Partial pressure of CO2 and CO2 mols
Calcite 0.0 0.185 dissolve only 

GAS_PHASE
    -fixed_pressure
    -pressure 100
    CO2(g) 1e-20
    H2S(g) 1e-20
    Hdg(g) 1e-20
    H2O(g) 1e-20

Use solution 1
REACTION 1
    Hdg 1
    0.06 in 10 steps
 

RUN_CELLS
-cell 1



SELECTED_OUTPUT
-reset false
-time true
-file saida.txt #Output file location goes here
-simulation true
-reaction true
-totals Na Cl Ca K C S Mg H H(0)
-equilibrium_phases Calcite
-temperature true
-gases CO2(g) CH4(g) H2S(g) Hdg(g)
-water
-charge_balance true
-ionic_strength true



USER_PUNCH
-headings Delta_Qtz Delta_Dol Delta_Cal Delta_Ill Delta_K-Feld Delta_Kao Delta_Chl Delta_Alb Si Ca Mg K Na Porosity
10 PUNCH KIN_DELTA("Quartz")
20 PUNCH KIN_DELTA("Dolomite")
30 PUNCH EQUI_DELTA("Calcite")
40 PUNCH KIN_DELTA("Illite")
50 PUNCH KIN_DELTA("K-Feldspar")
60 PUNCH KIN_DELTA("Kaolinite")
70 PUNCH KIN_DELTA("Chlorite14A")
80 PUNCH KIN_DELTA("Albite")
90 PUNCH TOT("Si")*1000*28.0843
100 PUNCH TOT("Ca")*1000*40.08
110 PUNCH TOT("Mg")*1000*24.312
120 PUNCH TOT("K")*1000*39.102
130 PUNCH TOT("Na")*1000*22.9898
REM Calculate Porosity
140 qtz = KIN("Quartz")
150 dol = KIN("Dolomite")
160 cal = EQUI("Calcite")
170 ill = KIN("Illite")
180 kao = KIN("Kaolinite")
190 fel = KIN("K-Feldspar")
200 alb = KIN("Albite")
210 qtz_V = (qtz * 22.67)
220 dol_V = (dol * 64.5)
230 cal_V = (cal * 36.9)
240 ill_V = (ill * 141.48)
250 kao_V = (kao * 99.35)
260 fel_V = (fel * 108.15)
270 alb_V = (alb * 101.31)
280 Vtot = 1000 + qtz_V + dol_V + cal_V + ill_V + kao_V + fel_V + alb_V
290 new_por = 1-(qtz_V + dol_V + cal_V + ill_V + kao_V + fel_V +alb_V)/Vtot
300 PUNCH new_por
-end

END

[dlparkhurst]

Sorry, I don't answer the kind of question you are asking. You will have to look at the output and plot results to make sure you understand the calculation and that it is doing what you want.

As for the reaction between Hdg and H2, that determines how fast reactive electron acceptors are added to your reactions. Ultimately, S(6) and C(4) will be consumed if enough H2 is added to the system.

[amaury22]

Dear Parkhurst

I also apologize as I made a very general question, but I actually had a specific doubt which was whether it is possible to use RUN_CELLS
-cell 1
and Kinetic at the same time, as I thought Run Cell was for transport, but it only has one cell, so maybe I was right. I repeat, my doubt is purely technical regarding Phreeqc

[dlparkhurst]

RUN_CELLS does no transport. It simply runs one or more cell batch reactions. For every cell listed in -cells, RUN_CELLS will define a reaction with the SOLUTION or MIX plus all of the reactants of the same number, including KINETICS,  and run the reaction. KINETICS will be run for the time specified by -time_step. After the reaction is calculated, the state of the SOLUTION and all other reactants will be SAVEd with the same user number.

[amaury22]

Dear Parkhurst

But I enter with the doubt that you are executing the actualization of the kinetics within the reaction, or if the kinetics is only considering that part initiated within the kinetics, there is no competition between the two?. Until now everything is working fine, I'm just trying to understand it well, so the calcite doesn't change,
just decrease Ca, but Calcite remains static

Thanks in advance Amaury

[dlparkhurst]

I am not completely sure what your question is.

You would normally not include Calcite as both EQUILIBRIUM_PHASES and KINETICS. If your time scale is long relative to the time it takes kinetic dissolution to reach equilibrium, use EQUILIBRIUM_PHASES; for short times, use KINETICS.

If you are talking about competition for H2, as you add H2, the most thermodynamically favorable electron acceptor will be used first, sulfate in your case. As more H2 is added, sulfate will be nearly completely reduced, and carbonate will begin to be reduced to methane.