MINERAL CHANGE WITH RESPECT TO TIME

[kennedyantwi1]

Good afternoon dear members,
I'm having challenges with simulating mineral changes in a caprock formation with respect to time. The code gives the following error messages:

ERROR: Elements in species have not been tabulated, H4SiO4.
ERROR: Reaction for species has not been defined, H4SiO4.
ERROR: Calculations terminating due to input errors.
Stopping.

I have enclosed the input script as below for your attention. Please, any guidelines, corrections and suggestions would be appreciated.
Thanks.
kennedyantwi01


# MODEL VALIDATION USING CO2-CAPROCK EXPERIMENTAL DATA
DATABASE c:\phreeqc\database\llnl.dat

SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
     log_k 16.681
     delta_h -5.738 kcal
    -analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9
    -dw 1.92e-9
    -Vm 7.29 0.92 2.07 -1.23 -1.60 # ref. 1 + McBride et al. 2015, JCED 60, 171
   
2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T
     log_k -1.8
    -analytical_expression 8.68 -0.0103 -2190
    -Vm 14.58 1.84 4.14 -2.46 -3.20
END

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
     
Chlorite(14A)
      Mg5Al2Si3O10(OH)8 + 16H+ = 5Mg+2 + 2Al+3 + 3H4SiO4 + 6H2O
      -log_k 68.38
      -delta_h -151.494 kcal
END
SOLUTION 1 # cap rock
-temperature 40.00
-pH 8.1
-pe 4.0
-redox pe
-units mg/l
Ca 142.9
Mg 55.2
Na 3777.0
Cl 4826.2
S(6) 785.9 # SO4--
C(4) 1410.9 # HCO3-
-density 1 calc
water 1 # kg
pressure 75.0 atm

EQUILIBRIUM_PHASES 1
CO2(g) 2.10 0.08345 # Partial pressure of CO2 and CO2 mols

EQUILIBRIUM_PHASES 1
Gypsum 0.0 0.0
K-feldspar 0.0 0.0 # microline
Siderite 0.0 0.0
Kaolinite 0.0 0.0
 
Save Solution 1
Save Equilibrium_Phases 1
END

RATES
Quartz
 -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_qtz = SI("Quartz")
30 if (M <= 0 and si_qtz < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_qtz > 0) then SA = 1e-05 #nucleation
60 k_acid = 0
70 k_neut = 10^(-13.99)*EXP(-87.60e+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-(10^si_qtz))
190 moles = r * TIME
200 SAVE moles
 -end
 
 Muscovite
 -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_musc = SI("Muscovite")
30 if (M <= 0 and si_musc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_musc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-11.85)*EXP(-22.0e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.370)
70 k_neut = 10^(-13.55)*EXP(-22.0e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-14.55)*EXP(-22.0e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.220)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_musc))
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.317))
70 k_neut = 10^(-12.56)*EXP(-65.0e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-66.50e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.471))
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 
   
 Chlorite-14A
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_hema = SI("Chlorite")
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 = 10^(-11.11)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.50)
70 k_neut = 10^(-12.52)*EXP(-88.00e+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-(10^si_chlo))
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   
   
 Dolomite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_dolo = SI("Dolomite")
30 if (M <= 0 and si_dolo < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_dolo > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_dolo))
190 moles = r * TIME
200 SAVE moles
    -end   
   
   
 KINETICS 1
Quartz
    -formula SiO2 1
    -m 60.08
    -m0 60.08
    -parms 1.36 # RSA in m2/mol
    -tol 1e-08
   
Muscovite
    -formula KAl3Si3O10(OH)2 1
    -m 398.71
    -m0 398.71
    -parms 0.00172
    -tol 1e-08

Albite
    -formula NaAlSi3O8 1
    -m 263.02
    -m0 263.02
    -parms 7.603e-5
    -tol 1e-08
 
Chlorite (14A) 
    -formula Mg5Al2Si3O10(OH)8 1 
    -m 341.76
    -m0 341.76
    -parms 0.002165
    -tol 1e-08   
   
Calcite
    -formula CaCO3 1
    -m 100.09
    -m0 100.09
    -parms 0.001499
    -tol 1e-08
   
Dolomite
     -formula CaMg(CO3)2 1
     -m 184.4
     -m0 184.4
     -parms 0.000976
     -tol 1e-08 
     -steps 3 30 300 3000 30000 300000 3000000 30000000 300000000 3000000000 30000000000 300000000000
    -cvode true
     
KNOBS
    -convergence_tolerance 1e-10
   
    INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings Time Quartz Muscovite Albite Chlorite Calcite Dolomite
    -axis_titles "Log10 Time" "Mineral (mol)"
    -initial_solutions false
    -connect_simulations true
    -plot_concentration_vs x
    -start
10 GRAPH_X log10(total_time)
20 GRAPH_Y kin("Quartz"), kin("Muscovite"), kin("Albite"), kin("Chlorite (14A)"), kin("Calcite"), kin("Dolomite")
    -end
   
-active true

SELECTED_OUTPUT 1
    -file POROSITY_18.05.2023_001.xls
    -reset false
    -time true
    -step true
    -ph true
    -pe true
    -totals C(4)
    -molalities K+ Na+ Ca+2 Mg+2 HCO3- SO4-2 Cl- CO3-2 CO2 H+ SiO2
    -kinetic_reactants Quartz Muscovite Albite Chlorite (14A) Calcite Dolomite
END

[dlparkhurst]

Please use the # button at the top of the Forum input screen to include code (and try to make your code as simple as possible).

H4SiO4 is not a species in the llnl.dat database. SiO2 is the master species. The following rewrites your equation in a way that will work with llnl.dat.

Code: [Select]

Chlorite(14A)
#      Mg5Al2Si3O10(OH)8 + 16H+ = 5Mg+2 + 2Al+3 + 3H4SiO4 + 6H2O
      Mg5Al2Si3O10(OH)8 + 16H+ = 5Mg+2 + 2Al+3 + 3SiO2 + 12H2O
      -log_k 68.38
      -delta_h -151.494 kcal


[kennedyantwi1]

Dear Dr. Parkhurst,
Thanks for you swift response. I applied your corrections and there were no error messages.
I am however not able to obtain the following output files:

1. Graphical Output of mineral changes
2. Excel output (xls) of the reaction minerals, brine composition, etc.

Please, I would appreciate your further assistance. The corrected code is attached for your attention.
Thanks in advance.
kennedyantwi01

Code: [Select]

# MODEL VALIDATION USING CO2-CAPROCK EXPERIMENTAL DATA
DATABASE c:\phreeqc\database\llnl.dat

SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
     log_k 16.681
     delta_h -5.738 kcal
    -analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9
    -dw 1.92e-9
    -Vm 7.29 0.92 2.07 -1.23 -1.60 # ref. 1 + McBride et al. 2015, JCED 60, 171
   
2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T
     log_k -1.8
    -analytical_expression 8.68 -0.0103 -2190
    -Vm 14.58 1.84 4.14 -2.46 -3.20
END

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
     
Chlorite(14A)
      Mg5Al2Si3O10(OH)8 + 16H+ = 5Mg+2 + 2Al+3 + 3SiO2 + 12H2O
      -log_k 68.38
      -delta_h -151.494 kcal
END
SOLUTION 1 # cap rock
-temperature 40.00
-pH 8.1
-pe 4.0
-redox pe
-units mg/l
Ca 142.9
Mg 55.2
Na 3777.0
Cl 4826.2
S(6) 785.9 # SO4--
C(4) 1410.9 # HCO3-
-density 1 calc
water 1 # kg
pressure 75.0 atm

EQUILIBRIUM_PHASES 1
CO2(g) 2.10 0.08345 # Partial pressure of CO2 and CO2 mols

EQUILIBRIUM_PHASES 1
Gypsum 0.0 0.0
K-feldspar 0.0 0.0 # microline
Siderite 0.0 0.0
Kaolinite 0.0 0.0
 
Save Solution 1
Save Equilibrium_Phases 1
END

RATES
Quartz
 -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_qtz = SI("Quartz")
30 if (M <= 0 and si_qtz < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_qtz > 0) then SA = 1e-05 #nucleation
60 k_acid = 0
70 k_neut = 10^(-13.99)*EXP(-87.60e+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-(10^si_qtz))
190 moles = r * TIME
200 SAVE moles
 -end
 
 Muscovite
 -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_musc = SI("Muscovite")
30 if (M <= 0 and si_musc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_musc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-11.85)*EXP(-22.0e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.370)
70 k_neut = 10^(-13.55)*EXP(-22.0e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-14.55)*EXP(-22.0e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.220)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_musc))
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.317))
70 k_neut = 10^(-12.56)*EXP(-65.0e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-66.50e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.471))
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 
   
 Chlorite-14A
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_hema = SI("Chlorite")
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 = 10^(-11.11)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.50)
70 k_neut = 10^(-12.52)*EXP(-88.00e+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-(10^si_chlo))
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   
   
 Dolomite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_dolo = SI("Dolomite")
30 if (M <= 0 and si_dolo < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_dolo > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_dolo))
190 moles = r * TIME
200 SAVE moles
    -end   
   
   
 KINETICS 1
Quartz
    -formula SiO2 1
    -m 60.08
    -m0 60.08
    -parms 1.36 # RSA in m2/mol
    -tol 1e-08
   
Muscovite
    -formula KAl3Si3O10(OH)2 1
    -m 398.71
    -m0 398.71
    -parms 0.00172
    -tol 1e-08

Albite
    -formula NaAlSi3O8 1
    -m 263.02
    -m0 263.02
    -parms 7.603e-5
    -tol 1e-08
 
Chlorite (14A) 
    -formula Mg5Al2Si3O10(OH)8 1 
    -m 341.76
    -m0 341.76
    -parms 0.002165
    -tol 1e-08   
   
Calcite
    -formula CaCO3 1
    -m 100.09
    -m0 100.09
    -parms 0.001499
    -tol 1e-08
   
Dolomite
     -formula CaMg(CO3)2 1
     -m 184.4
     -m0 184.4
     -parms 0.000976
     -tol 1e-08 
     -steps 3 30 300 3000 30000 300000 3000000 30000000 300000000 3000000000 30000000000 300000000000
    -cvode true
     
KNOBS
    -convergence_tolerance 1e-10
   
    INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings Time Quartz Muscovite Albite Chlorite Calcite Dolomite Gypsum K-feldspar Siderite   Kaolinite
    -axis_titles "Log10 Time" "Mineral (mol)"
    -initial_solutions false
    -connect_simulations true
    -plot_concentration_vs x
    -start
10 GRAPH_X log10(total_time)
20 GRAPH_Y kin("Quartz"), kin("Muscovite"), kin("Albite"), kin("Chlorite (14A)"), kin("Calcite"), kin("Dolomite"), kin("Gypsum"), kin("K-feldspar"), kin("Siderite"), kin("Kaolinite")
    -end 
-active true

PRINT
-reset false
-eh true
-equilibrium_phases true
-gas_phase true
-other true
-saturation_indices true
-species true
-totals true
-selected_output true

SELECTED_OUTPUT 1
    -file POROSITY_18.05.2023_001.xls
    -reset false
    -time true
    -step true
    -ph true
    -pe true
    -totals C(4)
    -molalities K+ Na+ Ca+2 Mg+2 HCO3- SO4-2 Cl- CO3-2 CO2 H+ Al+3 SiO2
    -kinetic_reactants Quartz Muscovite Albite Chlorite (14A) Calcite Dolomite
END


[dlparkhurst]

1. The Basic function KIN_DELTA give the change in the kinetic reaction (KIN) over the last integration period. KIN_TIME gives the length of time for the last integration, in seconds. KIN_DELTA("name") / KIN_TIME is the average rate in mol/sec for the last time period.

2. SELECTED_OUTPUT n allows you to write many quantities to a tab delimited file that can be read by Excel. If you want more control, you can use USER_PUNCH n in addition to SELECTED_OUTPUT n, to write Basic scripts with PUNCH statements to add columns to the selected-output file. You can have multiple SELECTED_OUTPUT/USER_PUNCH definitions, each identified by a unique user number (n).

[kennedyantwi1]

Thanks Dr. Parkhurst.