Dissolution of calcite during hydrogen storage in the presence of CO2 and ions

[MMMM]

Hello everyone. I am working on simulating hydrogen storage and I want to investigate the effect of different gas and ions on the solubility of hydrogen and the solubility of minerals. I have included the code I wrote below. The problem I have is that this code is very sensitive to gas values and runs just with specific values of gas. When it does run, the results fluctuate greatly over time. For example, after 10 days, all the calcite is dissolved, but after 20 days, only 0.01 of it is dissolved. I would appreciate any guidance on how to solve this problem.

Code: [Select]

[#DATA phreeqc.dat

SOLUTION 1
    pressure 1        # atm  Pressure
    temp 20               # degree Celsius
    -water 1
    units mol/kgw          # ppm (Resrvoir Salinity )
      Mg 1
      Cl 2
END

  GAS_PHASE 1
      -fixed_pressure
      -pressure 98.69  #100 bar
      CO2(g)  98.69
      H2S(g) 1e-20
      H2(g)  1e-20
      Hdg(g)  1e-20
      H2O(g) 1e-20
      N2(g) 1e-20
   
 END
 
 
 REACTION_PRESSURE 1
 1 1000 in 20


 RATES 1
 Calcite          #Phreeqc.DAT
 -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

 KINETICS 1
 Calcite        #Phreeqc.DAT
      -M0  10              # initial moles   
      -tol 1e-12
      -parms 1.67e5   0.6  # cm^2/mol calcite, exp factor
 
-steps 0 1 2 3 5 10 20 30 day
  -step_divide 200
  -runge_kutta 6
  -cvode true
  -bad_step_max 1000
  -cvode_order 5
  -cvode_steps 200
 


 reaction 2
 H2(g) 1
 100
 USE solution  1
 USE REACTION 2
 USE GAS_PHASE 1
 USE KINETICS 1
 
REACTION_PRESSURE 2
98.69
REACTION_TEMPERATURE 2
100
 

 USER_GRAPH 1
-headings  Temp
-axis_titles "Time (day)" "Amount of dissolved Calcite, mol"
-axis_scale x_axis auto auto auto auto
-axis_scale y_axis auto auto auto auto
initial_solutions true
-start
10 x = TOTAL_TIME/(86400)
20 GRAPH_X x
30 graph_y (-KIN_DELTA("Calcite"))
-end
   

USER_GRAPH 2
-headings Time  H2  CO2
-axis_titles "Time (day)" "H2 Solubility in Water, mol/kgw" "CO2 Solubility in Water, mol/kgw"
-axis_scale x_axis auto auto auto aut
-axis_scale y_axis auto auto auto aut
initial_solutions true
-start
10 x = TOTAL_TIME/86400
20 GRAPH_X x
30 GRAPH_Y   mol("H2")
40 GRAPH_SY  mol("CO2")
-end/code]

[dlparkhurst]

I don't really know where to start.

I guess the most important point is that the introduction of a large amount of H2(g) will lead to the complete reduction of C(4) to C(-4) (methane) in the solution. Whatever KINETICS/RATES definition is used, any dissolution of calcite will result in the formation of methane.

Code: [Select]

CaCO3 + 4H2 = Ca+2 + 2OH- + CH4 + H2O

The question really is how fast H2 reacts with C(4). You can model this by introducing the inert form Hdg(g), and then kinetically converting it to H2(g), which would react instantly. Thus, the Hdg to H2 kinetic reaction would determine the rate of C(4) reduction.

There are many posts related to hydrogen storage. Try searching for Hdg in the forum.

Here is a simplified script that demonstrates the problem with your simulation. Look at the output carefully.

Code: [Select]

SOLUTION 1
    pressure 1        # atm  Pressure
    temp 20               # degree Celsius
    -water 1
    units mol/kgw          # ppm (Resrvoir Salinity )
      Mg 1
      Cl 2

EQUILIBRIUM_PHASES 1
Calcite 0 1

GAS_PHASE 1
H2(g) 0

REACTION 1
 H2(g) 1
 100
END


[MMMM]

Thank you for kindly providing guidance.