help with initial condition for H2 reduction

[amaury22]

 I am relatively new to Phreeqc, and to learn, I gathered information from various sources and the manual to create my code. I wanted to study how H2 reacts in a medium with other minerals and salts. Although I don't encounter any errors, the output doesn't make much sense because all the values in the gas phase are zero. I believe it may have to do with the temperature and initial conditions, but I'm not sure how to proceed.

Code: [Select]

TITLE Mineral with H2

PHASES

H2(g)
        H2 = H2
        log_k -3.19
        delta_h 0.0

H2O(g)
        H2O = H2O
        log_k -4.27
        delta_h -9.82

CO2(g)
        CO2 = CO2
        log_k -2.87
        delta_h -3.94

H2S(g)
        H2S = H2S
        log_k 0.63
        delta_h -4.96

SOLUTION 1
    units mg/l
    temp 55
    pH 6.66
    density 1

    Na 104140 charge
    Ca 2365
    K 398
    Mg 836
    Sr 58
    Fe 2.5
    Al 5.8
    Cl 155610

    -water 1 # Defining water explicitly

EQUILIBRIUM_PHASES 1
    Quartz 0 51.06
    K-Feldspar 0 2.895         
    Hematite 0 3.537         
    Pyrite 0 3.826             
    Dolomite 0 1.457           
    Anhydrite 0 1.063         
    Calcite 0 1.375           
    Halite 0 0.961             
    Illite 0 0.358             
    Chlorite(14A) 0 0.0495     
    H2(g) 0 0       
    H2O(g) 0 0   


SAVE Solution 1 "solution.txt"
SAVE Equilibrium_Phases 1 "equilibrium_phases.txt"

USE Solution 1
USE Equilibrium_Phases 1


GAS_PHASE
    -fixed_pressure
    -pressure 81
    CO2(g) 1e-20
    H2S(g) 1e-20
    H2(g) 1e-20
    H2O(g) 1e-20

REACTION 1
    H2 1
    0 in 50 steps

SELECTED_OUTPUT
    -reset false
    -file graph_output.txt
    -simulation true
    -reaction true
    -totals Na Ca K Mg Sr Fe Al Cl H O
    -equilibrium_phases Quartz K-Feldspar Hematite Pyrite Dolomite Anhydrite Calcite Halite Illite Chlorite(14A)
    -si CO2(g) CH4(g) H2(g) H2O(g)
    -gas CO2(g) CH4(g) H2(g) H2O(g)


KNOBS
    -convergence_tolerance 1e-12
    -diagonal_scale 1e-8
    -iterations 1000

END

[dlparkhurst]

Don't define the gases with PHASES the way that you have. Use phreeqc.dat or pitzer.dat, which have Peng-Robinson definitions.

I would start by using Hdg(g) in place of H2(g), which will consider an unreactive H2 for initial calculations. You can see how much Hdg is needed to form a gas phase, and the solubility of Hdg (H2(g)) under your conditions.

I think that H2 will not react with your minerals very quickly, nor do I know whether H2 will reduce C(4) and S(6) at all, so you may want to consider a KINETIC reaction that converts inert Hdg to reactive H2. There are other posts on the forum that consider this approach.

[amaury22]

Dear Dlparkhurst

I do all that you suggest. I see the forum and incorporate the recumendation. However,
the problem of Hdg(g) zero in the outpu file still

Code: [Select]

DATA phreeqc.dat

TITLE H2-Gas_Mineral_Reaction_in_Mondarra_Gas_Reservoir

Selected_Output
   -file calciteout.txt
   -temperature true
   -totals C(4) C(-4) H(0) S(6) S(-2) Na Cl   # Certifique-se de incluir Hdg(g) aqui
   -kinetics Calcite
   -gases CO2(g) CH4(g) H2S(g) Hdg(g)
   -water
   -charge_balance true
   -ionic_strength true


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

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   
END

#================
# Mondarra Sample
#================
Solution 1
    -temperature           117.85
    -density                0.949
    -pressure               88.8231
    -pH                     10 charge
    units ppm
    Na                     47520
    Cl                     86500
    Ca                     5840
    K                      1152
    C                      596.7   as HCO3
    S(6)                   61.5 as SO4
    Mg                     771
    -water                 0.54
REACTION_TEMPERATURE 1
   117.85
REACTION_PRESSURE 1
   88.8231
   
Save Solution 1
END

Use Solution 1
KINETICS 1
Calcite
     -M0  6.5234893                    # initial moles   
     -tol 1e-8
     -parms 10 0.3

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

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

INCREMENTAL_REACTIONS true
Save Solution 2

END

the output is
  sim       reaction           temp             mu       mass_H2O         charge             Na             Ca              K             Mg             Sr             Fe             Al             Cl              H              O           C(4)          C(-4)           H(0)           S(6)          S(-2)             Na             Cl         Quartz       d_Quartz     K-Feldspar   d_K-Feldspar       Hematite     d_Hematite         Pyrite       d_Pyrite       Dolomite     d_Dolomite      Anhydrite    d_Anhydrite        Calcite      d_Calcite         Halite       d_Halite         Illite       d_Illite   Chlorite(14A)   d_Chlorite(14A)       pressure      total mol         volume       g_CO2(g)       g_CH4(g)       g_H2S(g)       g_Hdg(g)      k_Calcite     dk_Calcite   
           2            -99         55.000        6.55459       0.997597    2.71065e-11     6.1380e+00     6.8267e-02     1.8194e-01     2.2543e-02     9.0084e-04     1.7628e-06     2.7093e-07     6.4077e+00     9.6940e-12     0.0000e+00     1.1920e-03     9.4172e-14     9.6940e-12     4.7282e-02     1.7706e-09     6.1380e+00     6.4077e+00     5.1382e+01     3.2223e-01     2.6771e+00    -2.1793e-01     3.5369e+00    -8.1488e-05     3.8262e+00     2.2199e-04     1.3962e+00    -6.0807e-02     1.0154e+00    -4.7612e-02     1.4704e+00     1.2042e-01     5.2755e-01    -4.3345e-01     4.4174e-01     8.3736e-02     6.2313e-02     1.2813e-02     8.1000e+01     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     6.5235e+00    -0.0000e+00   
           2            -99         55.000        6.55459       0.997622    2.71084e-11     6.1380e+00     6.8266e-02     1.8194e-01     2.2543e-02     9.0082e-04     1.7628e-06     2.7093e-07     6.4077e+00     9.6940e-12     0.0000e+00     1.1920e-03     9.4172e-14     9.6940e-12     4.7282e-02     1.7706e-09     6.1380e+00     6.4077e+00     5.1382e+01    -8.7792e-05     2.6771e+00     6.2135e-05     3.5368e+00    -1.0961e-04     3.8264e+00     2.1922e-04     1.3958e+00    -4.4004e-04     1.0149e+00    -4.4048e-04     1.4713e+00     8.8005e-04     5.2740e-01    -1.5749e-04     4.4163e-01    -1.0822e-04     6.2407e-02     9.3383e-05     8.1000e+01     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     6.5235e+00    -0.0000e+00   
           2            -99         55.000        6.55459       0.997647    2.71061e-11     6.1380e+00     6.8265e-02     1.8194e-01     2.2542e-02     9.0079e-04     1.7627e-06     2.7093e-07     6.4077e+00     9.6940e-12     0.0000e+00     1.1920e-03     9.4172e-14     9.6940e-12     4.7283e-02     1.7706e-09     6.1380e+00     6.4077e+00     5.1382e+01    -8.7793e-05     2.6772e+00     6.2135e-05     3.5367e+00    -1.0961e-04     3.8267e+00     2.1922e-04     1.3953e+00    -4.4004e-04     1.0145e+00    -4.4048e-04     1.4722e+00     8.8005e-04     5.2724e-01    -1.5749e-04     4.4152e-01    -1.0822e-04     6.2500e-02     9.3383e-05     8.1000e+01     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     6.5235e+00    -0.0000e+00   
           2            -99         55.000        6.55458       0.997671    2.71018e-11     6.1380e+00     6.8264e-02     1.8194e-01     2.2542e-02     9.0077e-04     1.7627e-06     2.7093e-07     6.4077e+00     9.6940e-12     0.0000e+00     1.1920e-03     9.4172e-14     9.6940e-12     4.7284e-02     1.7707e-09     6.1380e+00     6.4077e+00     5.1382e+01    -8.7793e-05     2.6773e+00     6.2135e-05     3.5366e+00    -1.0961e-04     3.8269e+00     2.1922e-04     1.3949e+00    -4.4004e-04     1.0141e+00    -4.4048e-04     1.4731e+00     8.8005e-04     5.2708e-01    -1.5749e-04     4.4141e-01    -1.0822e-04     6.2594e-02     9.3383e-05     8.1000e+01     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     6.5235e+00    -0.0000e+00   
           2            -99         55.000        6.55458       0.997696    2.71013e-11     6.1381e+00     6.8262e-02     1.8194e-01     2.2541e-02     9.0075e-04     1.7627e-06     2.7094e-07     6.4077e+00     9.6940e-12     0.0000e+00     1.1920e-03     9.4172e-14     9.6940e-12     4.7285e-02     1.7707e-09     6.1381e+00     6.4077e+00     5.1382e+01    -8.7793e-05     2.6773e+00     6.2135e-05     3.5365e+00    -1.0961e-04     3.8271e+00     2.1922e-04     1.3944e+00    -4.4004e-04     1.0136e+00    -4.4048e-04     1.4739e+00     8.8005e-04     5.2692e-01    -1.5749e-04     4.4130e-01    -1.0822e-04     6.2687e-02     9.3383e-05     8.1000e+01     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     0.0000e+00     6.5235e+00    -0.0000e+00   

[dlparkhurst]

Here is something to get you started. It adds Hdg to your solution at the T and P you specify. After adding about 0.035 mol of Hdg to the system, a gas phase forms.

Assuming that H2 does not react, this calculation describes your system. If H2 does react, you can use kinetics -formula Hdg -1 H2 +1 to begin to react with the minerals. I'm not convinced that H2 will react with the minerals in this way, but it is your calculation. I would also assume that calcite is in equilibrium rather than treating it kinetically.

Code: [Select]

Solution 1
    -temperature           117.85
    -density                0.949
    -pressure               88.8231
    -pH                     10 charge
    units ppm
    Na                     47520
    Cl                     86500
    Ca                     5840
    K                      1152
    C                      596.7   as HCO3
    S(6)                   61.5 as SO4
    Mg                     771
    #-water                 0.54
REACTION_TEMPERATURE 1
   117.85
REACTION_PRESSURE 1
   88.8231
GAS_PHASE
Hdg(g) 0
REACTION
Hdg 1
.06 mole in 30
USER_GRAPH 1
    -headings               rxn TOT(Hdg) GAS(Hdg)
    -axis_titles            "Hdg added, moles" "Dissolved Hdg, mol/kgw" "Gas Hdg, moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN
20 GRAPH_Y TOT("Hdg")
30 GRAPH_SY GAS("Hdg(g)")
  -end
    -active                 true
END


[amaury22]

Dear Dlparkhurst

Thanks a lot. This exemple is functioning well with your recomendation. Now, we do the following step to add If H2 does react, you can use kinetics -formula Hdg -1 H2 +1 to begin to react with the minerals

Code: [Select]

DATA phreeqc.dat

TITLE Mineral and Hdg

SOLUTION 1
    pressure 81
    temp 55
    pH 6.66
    -density 1 calc
    -units mg/l
    Na 104140 charge
    Ca 2365
    K 398
    Mg 836
    Sr 58
    Fe 2.5
    Al 5.8
    Cl 155610


EQUILIBRIUM_PHASES 1
    Quartz 0 51.06
    K-Feldspar 0 2.895         
    Hematite 0 3.537         
    Pyrite 0 3.826             
    Dolomite 0 1.457           
    Anhydrite 0 1.063         
    Calcite 0 1.375           
    Halite 0 0.961             
    Illite 0 0.358             
    Chlorite(14A) 0 0.0495     
   

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

REACTION 1
    Hdg 1
    0.03 in 30 steps

INCREMENTAL_REACTIONS true

SELECTED_OUTPUT
    -reset false
    -file graph_output.txt
    -simulation true
    -reaction true
    -totals Na Cl Ca K C S Mg
    -equilibrium_phases Quartz K-Feldspar Hematite Pyrite Dolomite Anhydrite Calcite Halite Illite Chlorite(14A)
    -temperature true
    -gases CO2(g) CH4(g) H2S(g) Hdg(g)
    -water
    -charge_balance true
    -ionic_strength true

KNOBS
    -convergence_tolerance 1e-14
    -diagonal_scale 1e-8
    -iterations 31000

END

[amaury22]

Caro Parkhurst

I have done all the recommendations you gave me and everything seems to be working fine. Except for a certain warning that I don't know if it's important, but it seems important because it starts giving non-zero values after step 20.

What do you consider?

Thanks for everything

o warning is:
Simulation 1. Kinetic step 18.                                                 
WARNING: Maximum iterations exceeded, 31000

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying smaller step size, pe step size 10, 5 ...


Simulation 1. Kinetic step 19.                                                 
WARNING: Maximum iterations exceeded, 62000

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying reduced tolerance 1e-16 ...

several warning like this

Code: [Select]

DATA phreeqc.dat

TITLE Mineral and Hdg

SOLUTION 1
    pressure 81
    temp 55
    pH 6.66
    -density 1 calc
    -units mg/l
    Na 104140 charge
    Ca 2365
    K 398
    Mg 836
    Sr 58
    Fe 2.5
    Al 5.8
    Cl 155610

EQUILIBRIUM_PHASES 1
    Quartz 0 51.06
    K-Feldspar 0 2.895         
    Hematite 0 3.537         
    Pyrite 0 3.826             
    Dolomite 0 1.457           
    Anhydrite 0 1.063         
    Calcite 0 1.375           
    Halite 0 0.961             
    Illite 0 0.358             
    Chlorite(14A) 0 0.0495     

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

REACTION 1
    Hdg 1
    0.03 in 30 steps

KINETICS 1
    -formula Hdg -1 H2 +1
    -m0 1.0
    -steps 30
    -time_steps 1 day

INCREMENTAL_REACTIONS true

SELECTED_OUTPUT
    -reset false
    -file graph_output.txt
    -simulation true
    -reaction true
    -totals Na Cl Ca K C S Mg
    -equilibrium_phases Quartz K-Feldspar Hematite Pyrite Dolomite Anhydrite Calcite Halite Illite Chlorite(14A)
    -temperature true
    -gases CO2(g) CH4(g) H2S(g) Hdg(g)
    -water
    -charge_balance true
    -ionic_strength true

KNOBS
    -convergence_tolerance 1e-14
    -diagonal_scale 1e-8
    -iterations 31000

END

[dlparkhurst]

I'm not sure how it is doing anything. You need a RATES definition to go with your kinetic reaction. However, I suggest you simply add H2 in REACTION an look at how the minerals react.

Code: [Select]

DATA phreeqc.dat

TITLE Mineral and Hdg

SOLUTION 1
    pressure 81
    temp 55
    pH 6.66
    -density 1 calc
    -units mg/l
    Na 104140 charge
    Ca 2365
    K 398
    Mg 836
    Sr 58
    Fe 2.5
    Al 5.8
    Cl 155610

EQUILIBRIUM_PHASES 1
    Quartz 0 51.06
    K-Feldspar 0 2.895         
    Hematite 0 3.537         
    Pyrite 0 3.826             
    Dolomite 0 1.457           
    Anhydrite 0 1.063         
    Calcite 0 1.375           
    Halite 0 0.961             
    Illite 0 0.358             
    Chlorite(14A) 0 0.0495     
SAVE solution 2
SAVE equilibrium_phases 2
END
USE solution 2
USE equilibrium_phases 2
REACTION 1
    H2 1
    0.01
END

[amaury22]

Dear Parkhurst

I'm sorry, I sent the wrong code, look at the one I put now.
As I told you, I did what you said and it seems to be working fine.
It only has several warnings of the type
Simulation 1. Kinetic step 20.
WARNING: Maximum iterations exceeded, 62000

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying to scale pure_phase columns 1e-10 ...


Simulation 1. Kinetic step 20.
WARNING: Maximum iterations exceeded, 62000

WARNING: Numerical method failed with this set of convergence parameters.

WARNING: Trying to scale pure_phase columns and diagonal scale 1e-10 ...

Code: [Select]

DATA phreeqc.dat

TITLE Mineral and Hdg

SOLUTION 1
    pressure 81
    temp 55
    pH 6.66
    -density 1 calc
    -units mg/l
    Na 104140 charge
    Ca 2365
    K 398
    Mg 836
    Sr 58
    Fe 2.5
    Al 5.8
    Cl 155610

EQUILIBRIUM_PHASES 1
       
 Calcite 0 2.375           
     

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

REACTION 1
    Hdg 1
    0.03 in 30 steps


RATES
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

Use Solution 1
KINETICS 1

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

-bad_step_max 55000


INCREMENTAL_REACTIONS true

SELECTED_OUTPUT
    -reset false
    -file graph_output.txt
    -simulation true
    -reaction true
    -totals Na Cl Ca K C S Mg
    -equilibrium_phases Calcite
    -temperature true
    -gases CO2(g) CH4(g) H2S(g) Hdg(g)
    -water
    -charge_balance true
    -ionic_strength true

KNOBS
    -convergence_tolerance 1e-14
    -diagonal_scale 1e-8
    -iterations 31000

END

[dlparkhurst]

This is the KINETICS version of my last post, but I still think that the last script is more than enough for you to consider.

Code: [Select]

#DATA phreeqc.dat

TITLE Mineral and Hdg

SOLUTION 1
    pressure 81
    temp 55
    pH 6.66
    -density 1 calc
    -units mg/l
    Na 104140 charge
    Ca 2365
    K 398
    Mg 836
    Sr 58
    Fe 2.5
    Al 5.8
    Cl 155610
END
EQUILIBRIUM_PHASES 1
       
 Calcite 0 2.375           
     
END
GAS_PHASE 1
    -fixed_pressure
    -pressure 81
    CO2(g) 1e-20
    H2S(g) 1e-20
    Hdg(g) 1e-20
    H2O(g) 1e-20
END
REACTION 1
    Hdg 1
    0.03 #in 30 steps
END
RUN_CELLS
-cell 1
END

RATES
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 1

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
END
RUN_CELLS
-cell 1

END


[amaury22]

Dear Parkhurst

It code is functioning well without erros, but calcite not change in the kinetic with Hdg. It is happening that Hdg hydrogen does not react with Calcite, is that it?
And how could this effect be considered by adding kinetic in calcite as well.
Yet you said that this was not good?

REgards,
Amaury