Saturation Index

[dat]

Hi,

I am finding the saturation index of minerals with the interaction of CO2 and brine.
Initially, I am going to equilibrate my minerals with brine and the CO2 will then inject.

But the issue is I am not getting the values in the paper, in their model and experiment.

The main issue I identified is in my llnl database mineral Chlorite and Montmorillonite is not there. But they have considered that for the model.

And also I hereby attach my code. Could you please if this is ok? Really appreciate the support.

Code: [Select]

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 = (1e-14)*EXP(-87.7e+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_qtz))
    190 moles = r * TIME
    200 SAVE moles
     -end
Anorthite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anort = SI("Anorthite")
   30 if (M <= 0 and si_anort < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_anort > 0) then SA = 1e-05  #nucleation
   60 k_acid = (1e-08)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))                 
   70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(si_anort))*ACT("H+")^(1.411)
   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 = (1e-12)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(si_alb))*(ACT("H+")^(0.5))
    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))               
    70 k_neut = 0
    80 k_carb = 0
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(si_calc))*ACT("H+")^(0.5)
    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 = (3.98e-13)*EXP(-46.00e+03/8.314*(1.0/TK-1.0/373.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(si_ill))*ACT("H+")^(0.1)
    190 moles = r * TIME
    200 SAVE moles
    -end

Clinochlore-14A
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_chlo = SI("Clinochlore-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
K-feldspar
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_felds = SI("K-feldspar")
   30 if (M <= 0 and si_felds < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_felds > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = (8.71e-11)*EXP(-51.7e+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_felds))*ACT("H+")^0.5
   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
Pyrite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_py = SI("Pyrite")
   30 if (M <= 0 and si_py < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_py > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = (4.0e-11)*EXP(-62.76e+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_py))*ACT("H+")^0.5
   190 moles = r * TIME
   200 SAVE moles
-end
Goethite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_goe = SI("Goethite")
   30 if (M <= 0 and si_goe < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_goe > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = (1.15e-13)*EXP(-86e+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_goe))*ACT("H+")^0.5
   190 moles = r * TIME
   200 SAVE moles
-end

SOLUTION 1
    temp      40
    pH        4.41
    pe        2
    redox     pe
    units     ppm
    density   2.1 calc
    Cl        27335
    Na        17710
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Albite    0 0.18
    Anorthite 0 0.17
    Goethite  0 0.54
    Illite    0 0.47
    K-Feldspar 0 0.26
    Kaolinite 0 0.19
    Pyrite    0 0.07
    Quartz    0 4.26

SAVE SOLUTION 1
SAVE Equilibrium_Phases 1
END     

USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 40
    CO2(g)    100

INCREMENTAL_REACTIONS True


KINETICS 1
Quartz
    -formula  SiO2  1
    -m        4.26
    -m0       4.26
    -parms    2.16
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.18
    -m0       0.18
    -parms    9.5469
    -tol      1e-08
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        0.17
    -m0       0.17
    -parms    9.6188
    -tol      1e-08
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        0.47
    -m0       0.47
    -parms    179.5235
    -tol      1e-08
K-feldspar
    -formula  K-feldspar  1
    -m        0.26
    -m0       0.26
    -parms    10.37
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        0.07
    -m0       0.07
    -parms    4.32
    -tol      1e-08
Goethite
    -formula  FeOOH  1
    -m        0.54
    -m0       0.54
    -parms    1.958
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0.19
    -m0       0.19
    -parms    298.377
    -tol      1e-08

-steps       14 35 63 126 259 day
-step_divide 1
-runge_kutta 3
-bad_step_max 500



[dlparkhurst]

The databases distributed with PHREEQC have various formulations for chlorite (14A and 7A, referring to the layer spacing) and montmorillionite (Ca-Montmorillonite, Belle Fourche, Aberdeen and others. Look at phreeqc.dat, wateq4f.dat, sit.dat, and minteq.dat. Hopefully, the paper gives the formula and log Ks fro the mineral compositions that they used. If so, you can define them with PHASES in your input file.

I resist "checking" your script. I do not know your conceptual model so you will have to examine the output and decide if it fits your expectations.

[dat]

Hi,

Thank you. I tried adding the new phases.

Code: [Select]

SOLUTION_MASTER_SPECIES

Si H4SiO4 0 SiO2 28.0843

SOLUTION_SPECIES
H4SiO4 = H4SiO4
-dw 1.10e-9
-Vm  10.5  1.7  20  -2.7  0.1291 # supcrt + 2*H2O in a1

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

Montmorillonite-BCCa
Ca0.17Mg0.34Al1.66Si4O10(OH)2     = 0.170Ca+2     + 0.340Mg+2     + 1.660Al+3     - 6.000H+     + 4.000H4(SiO4)     - 4.000H2O
     log_k     4.200     
     delta_h -156.000    #kJ/mol   
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 = (1e-14)*EXP(-87.7e+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_qtz))
    190 moles = r * TIME
    200 SAVE moles
     -end
Anorthite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anort = SI("Anorthite")
   30 if (M <= 0 and si_anort < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_anort > 0) then SA = 1e-05  #nucleation
   60 k_acid = (1e-08)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))                 
   70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(si_anort))*ACT("H+")^(1.411)
   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 = (1e-12)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(si_alb))*(ACT("H+")^(0.5))
    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))               
    70 k_neut = 0
    80 k_carb = 0
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(si_calc))*ACT("H+")^(0.5)
    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 = (3.98e-13)*EXP(-46.00e+03/8.314*(1.0/TK-1.0/373.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(si_ill))*ACT("H+")^(0.1)
    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
K-feldspar
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_felds = SI("K-feldspar")
   30 if (M <= 0 and si_felds < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_felds > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = (8.71e-11)*EXP(-51.7e+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_felds))*ACT("H+")^0.5
   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
Pyrite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_py = SI("Pyrite")
   30 if (M <= 0 and si_py < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_py > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = (4.0e-11)*EXP(-62.76e+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_py))*ACT("H+")^0.5
   190 moles = r * TIME
   200 SAVE moles
-end
Goethite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_goe = SI("Goethite")
   30 if (M <= 0 and si_goe < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_goe > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = (1.15e-13)*EXP(-86e+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_goe))*ACT("H+")^0.5
   190 moles = r * TIME
   200 SAVE moles
-end
Montmorillonite-BCCa
-start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_mont = SI("Montmorillonite-BCCa")
   30 if (M <= 0 and si_mont < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_mont > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = (1.7e-13)*EXP(-35e+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_mont))*ACT("H+")^0.5
   190 moles = r * TIME
   200 SAVE moles
-end



SOLUTION 1
    temp      40
    pH        4.41
    pe        2
    redox     pe
    units     ppm
    density   2.1 calc
    Cl        27335
    Na        17710
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Albite    0 0.18
    Anorthite 0 0.17
    Goethite  0 0.54
    Illite    0 0.47
    K-Feldspar 0 0.26
    Kaolinite 0 0.19
    Pyrite    0 0.07
    Quartz    0 4.26

SAVE SOLUTION 1
SAVE Equilibrium_Phases 1
END     

USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1

GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 40
    CO2(g)    100

INCREMENTAL_REACTIONS True


KINETICS 1
Quartz
    -formula  SiO2  1
    -m        4.26
    -m0       4.26
    -parms    2.16
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.18
    -m0       0.18
    -parms    9.5469
    -tol      1e-08
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        0.17
    -m0       0.17
    -parms    9.6188
    -tol      1e-08
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        0.47
    -m0       0.47
    -parms    179.5235
    -tol      1e-08
K-feldspar
    -formula  K-feldspar  1
    -m        0.26
    -m0       0.26
    -parms    10.37
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        0.07
    -m0       0.07
    -parms    4.32
    -tol      1e-08
Goethite
    -formula  FeOOH  1
    -m        0.54
    -m0       0.54
    -parms    1.958
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0.19
    -m0       0.19
    -parms    298.377
    -tol      1e-08
Chlorite(14A)
    -formula  Mg5Al2Si3O10(OH)8  1
    -m        1
    -m0       1
    -parms    1
    -tol      1e-08
Montmorillonite-BCCa
    -formula  Ca0.17Mg0.34Al1.66Si4O10(OH)2  1
    -m        1
    -m0       1
    -parms    1
    -tol      1e-08
-steps       1209600 3024000 5443200 10886400 22377600
-step_divide 1
-runge_kutta 3
-bad_step_max 500


But I am getting some errors. Could you please assist me with this?

Code: [Select]

Reading data base.
------------------

LLNL_AQUEOUS_MODEL_PARAMETERS
NAMED_EXPRESSIONS
SOLUTION_MASTER_SPECIES
SOLUTION_SPECIES
PHASES
EXCHANGE_MASTER_SPECIES
EXCHANGE_SPECIES
SURFACE_MASTER_SPECIES
SURFACE_SPECIES
RATES
END
------------------------------------
Reading input data for simulation 1.
------------------------------------

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\llnl.dat
SOLUTION_MASTER_SPECIES
Si H4SiO4 0 SiO2 28.0843
SOLUTION_SPECIES
H4SiO4 = H4SiO4
dw 1.10e-9
vm  10.5  1.7  20  -2.7  0.1291 # supcrt + 2*H2O in a1
END
ERROR: Could not reduce equation to secondary master species, H2SiO4-2.
ERROR: Non-master species in secondary reaction, H2SiO4-2.
ERROR: Could not reduce equation to secondary master species, H4(H2SiO4)4-4.
ERROR: Non-master species in secondary reaction, H4(H2SiO4)4-4.
ERROR: Could not reduce equation to secondary master species, H6(H2SiO4)4-2.
ERROR: Non-master species in secondary reaction, H6(H2SiO4)4-2.
ERROR: Could not reduce equation to secondary master species, HSiO3-.
ERROR: Non-master species in secondary reaction, HSiO3-.
ERROR: Could not reduce equation to secondary master species, NaHSiO3.
ERROR: Non-master species in secondary reaction, NaHSiO3.
ERROR: Could not reduce equation to secondary master species, SiF6-2.
ERROR: Non-master species in secondary reaction, SiF6-2.
ERROR: Could not reduce equation to secondary master species, SiO2.
ERROR: Non-master species in secondary reaction, SiO2.
ERROR: Could not reduce equation to secondary master species, Afwillite.
ERROR: Could not reduce equation to secondary master species, Akermanite.
ERROR: Could not reduce equation to secondary master species, Alamosite.
ERROR: Could not reduce equation to secondary master species, Albite.
ERROR: Could not reduce equation to secondary master species, Albite_high.
ERROR: Could not reduce equation to secondary master species, Albite_low.
ERROR: Could not reduce equation to secondary master species, Amesite-14A.
ERROR: Could not reduce equation to secondary master species, Analcime.
ERROR: Could not reduce equation to secondary master species, Analcime-dehy.
ERROR: Could not reduce equation to secondary master species, Andalusite.
ERROR: Could not reduce equation to secondary master species, Andradite.
ERROR: Could not reduce equation to secondary master species, Annite.
ERROR: Could not reduce equation to secondary master species, Anorthite.
ERROR: Could not reduce equation to secondary master species, Anthophyllite.
ERROR: Could not reduce equation to secondary master species, Antigorite.
ERROR: Could not reduce equation to secondary master species, Ba2Si3O8.
ERROR: Could not reduce equation to secondary master species, Ba2SiO4.
ERROR: Could not reduce equation to secondary master species, BaSiF6.
ERROR: Could not reduce equation to secondary master species, Beidellite-Ca.
ERROR: Could not reduce equation to secondary master species, Beidellite-Cs.
ERROR: Could not reduce equation to secondary master species, Beidellite-H.
ERROR: Could not reduce equation to secondary master species, Beidellite-K.
ERROR: Could not reduce equation to secondary master species, Beidellite-Mg.
ERROR: Could not reduce equation to secondary master species, Beidellite-Na.
ERROR: Could not reduce equation to secondary master species, Boltwoodite.
ERROR: Could not reduce equation to secondary master species, Boltwoodite-Na.
ERROR: Could not reduce equation to secondary master species, Ca-Al_Pyroxene.
ERROR: Could not reduce equation to secondary master species, CdSiO3.
ERROR: Could not reduce equation to secondary master species, Celadonite.
ERROR: Could not reduce equation to secondary master species, Chalcedony.
ERROR: Could not reduce equation to secondary master species, Chamosite-7A.
ERROR: Could not reduce equation to secondary master species, Chrysocolla.
ERROR: Could not reduce equation to secondary master species, Chrysotile.
ERROR: Could not reduce equation to secondary master species, Clinochlore-14A.
ERROR: Could not reduce equation to secondary master species, Clinochlore-7A.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-Ca.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-Cs.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-dehy.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-dehy-Ca.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-dehy-Cs.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-dehy-K.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-dehy-Na.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-dehy-NH4.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-dehy-Sr.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-hy-Ca.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-hy-Cs.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-hy-K.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-hy-Na.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-hy-Sr.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-K.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-Na.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-NH4.
ERROR: Could not reduce equation to secondary master species, Clinoptilolite-Sr.
ERROR: Could not reduce equation to secondary master species, Clinozoisite.
ERROR: Could not reduce equation to secondary master species, Co2SiO4.
ERROR: Could not reduce equation to secondary master species, Coesite.
ERROR: Could not reduce equation to secondary master species, Coffinite.
ERROR: Could not reduce equation to secondary master species, Cordierite_anhyd.
ERROR: Could not reduce equation to secondary master species, Cordierite_hydr.
ERROR: Could not reduce equation to secondary master species, Cristobalite(alpha).
ERROR: Could not reduce equation to secondary master species, Cristobalite(beta).
ERROR: Could not reduce equation to secondary master species, Cronstedtite-7A.
ERROR: Could not reduce equation to secondary master species, Daphnite-14A.
ERROR: Could not reduce equation to secondary master species, Daphnite-7A.
ERROR: Could not reduce equation to secondary master species, Dicalcium_silicate.
ERROR: Could not reduce equation to secondary master species, Diopside.
ERROR: Could not reduce equation to secondary master species, Dioptase.
ERROR: Could not reduce equation to secondary master species, Enstatite.
ERROR: Could not reduce equation to secondary master species, Epidote.
ERROR: Could not reduce equation to secondary master species, Epidote-ord.
ERROR: Could not reduce equation to secondary master species, Eucryptite.
ERROR: Could not reduce equation to secondary master species, Fayalite.
ERROR: Could not reduce equation to secondary master species, Ferrosilite.
ERROR: Could not reduce equation to secondary master species, Forsterite.
ERROR: Could not reduce equation to secondary master species, Foshagite.
ERROR: Could not reduce equation to secondary master species, Gehlenite.
ERROR: Could not reduce equation to secondary master species, Gismondine.
ERROR: Could not reduce equation to secondary master species, Greenalite.
ERROR: Could not reduce equation to secondary master species, Grossular.
ERROR: Could not reduce equation to secondary master species, Gyrolite.
ERROR: Could not reduce equation to secondary master species, Haiweeite.
ERROR: Could not reduce equation to secondary master species, Hatrurite.
ERROR: Could not reduce equation to secondary master species, Hedenbergite.
ERROR: Could not reduce equation to secondary master species, Heulandite.
ERROR: Could not reduce equation to secondary master species, Hillebrandite.
ERROR: Could not reduce equation to secondary master species, Illite.
ERROR: Could not reduce equation to secondary master species, Jadeite.
ERROR: Could not reduce equation to secondary master species, K-Feldspar.
ERROR: Could not reduce equation to secondary master species, Kalsilite.
ERROR: Could not reduce equation to secondary master species, Kaolinite.
ERROR: Could not reduce equation to secondary master species, Kasolite.
ERROR: Could not reduce equation to secondary master species, Kyanite.
ERROR: Could not reduce equation to secondary master species, Larnite.
ERROR: Could not reduce equation to secondary master species, Laumontite.
ERROR: Could not reduce equation to secondary master species, Lawsonite.
ERROR: Could not reduce equation to secondary master species, Margarite.
ERROR: Could not reduce equation to secondary master species, Maximum_Microcline.
ERROR: Could not reduce equation to secondary master species, Merwinite.
ERROR: Could not reduce equation to secondary master species, Mesolite.
ERROR: Could not reduce equation to secondary master species, Minnesotaite.
ERROR: Could not reduce equation to secondary master species, Monticellite.
ERROR: Could not reduce equation to secondary master species, Montmor-Ca.
ERROR: Could not reduce equation to secondary master species, Montmor-Cs.
ERROR: Could not reduce equation to secondary master species, Montmor-K.
ERROR: Could not reduce equation to secondary master species, Montmor-Mg.
ERROR: Could not reduce equation to secondary master species, Montmor-Na.
ERROR: Could not reduce equation to secondary master species, Mordenite.
ERROR: Could not reduce equation to secondary master species, Mordenite-dehy.
ERROR: Could not reduce equation to secondary master species, Muscovite.
ERROR: Could not reduce equation to secondary master species, Na2SiO3.
ERROR: Could not reduce equation to secondary master species, Na4SiO4.
ERROR: Could not reduce equation to secondary master species, Na6Si2O7.
ERROR: Could not reduce equation to secondary master species, Natrolite.
ERROR: Could not reduce equation to secondary master species, Natrosilite.
ERROR: Could not reduce equation to secondary master species, Nepheline.
ERROR: Could not reduce equation to secondary master species, Ni2SiO4.
ERROR: Could not reduce equation to secondary master species, Nontronite-Ca.
ERROR: Could not reduce equation to secondary master species, Nontronite-Cs.
ERROR: Could not reduce equation to secondary master species, Nontronite-H.
ERROR: Could not reduce equation to secondary master species, Nontronite-K.
ERROR: Could not reduce equation to secondary master species, Nontronite-Mg.
ERROR: Could not reduce equation to secondary master species, Nontronite-Na.
ERROR: Could not reduce equation to secondary master species, Okenite.
ERROR: Could not reduce equation to secondary master species, Paragonite.
ERROR: Could not reduce equation to secondary master species, Pargasite.
ERROR: Could not reduce equation to secondary master species, Pb2SiO4.
ERROR: Could not reduce equation to secondary master species, Petalite.
ERROR: Could not reduce equation to secondary master species, Phlogopite.
ERROR: Could not reduce equation to secondary master species, Prehnite.
ERROR: Could not reduce equation to secondary master species, Pseudowollastonite.
ERROR: Could not reduce equation to secondary master species, Pyrophyllite.
ERROR: Could not reduce equation to secondary master species, Quartz.
ERROR: Could not reduce equation to secondary master species, Rankinite.
ERROR: Could not reduce equation to secondary master species, Rhodonite.
ERROR: Could not reduce equation to secondary master species, Ripidolite-14A.
ERROR: Could not reduce equation to secondary master species, Ripidolite-7A.
ERROR: Could not reduce equation to secondary master species, Sanbornite.
ERROR: Could not reduce equation to secondary master species, Sanidine_high.
ERROR: Could not reduce equation to secondary master species, Saponite-Ca.
ERROR: Could not reduce equation to secondary master species, Saponite-Cs.
ERROR: Could not reduce equation to secondary master species, Saponite-H.
ERROR: Could not reduce equation to secondary master species, Saponite-K.
ERROR: Could not reduce equation to secondary master species, Saponite-Mg.
ERROR: Could not reduce equation to secondary master species, Saponite-Na.
ERROR: Could not reduce equation to secondary master species, Scolecite.
ERROR: Could not reduce equation to secondary master species, Sepiolite.
ERROR: Could not reduce equation to secondary master species, Si.
ERROR: Could not reduce equation to secondary master species, Si(g).
ERROR: Could not reduce equation to secondary master species, SiF4(g).
ERROR: Could not reduce equation to secondary master species, Sillimanite.
ERROR: Could not reduce equation to secondary master species, SiO2(am).
ERROR: Could not reduce equation to secondary master species, Sklodowskite.
ERROR: Could not reduce equation to secondary master species, Smectite-high-Fe-Mg.
ERROR: Could not reduce equation to secondary master species, Smectite-low-Fe-Mg.
ERROR: Could not reduce equation to secondary master species, Soddyite.
ERROR: Could not reduce equation to secondary master species, Spodumene.
ERROR: Could not reduce equation to secondary master species, Sr2SiO4.
ERROR: Could not reduce equation to secondary master species, SrSiO3.
ERROR: Could not reduce equation to secondary master species, Stilbite.
ERROR: Could not reduce equation to secondary master species, Talc.
ERROR: Could not reduce equation to secondary master species, Tephroite.
ERROR: Could not reduce equation to secondary master species, Titanite.
ERROR: Could not reduce equation to secondary master species, Tobermorite-11A.
ERROR: Could not reduce equation to secondary master species, Tobermorite-14A.
ERROR: Could not reduce equation to secondary master species, Tobermorite-9A.
ERROR: Could not reduce equation to secondary master species, Tremolite.
ERROR: Could not reduce equation to secondary master species, Tridymite.
ERROR: Could not reduce equation to secondary master species, Uranophane.
ERROR: Could not reduce equation to secondary master species, Wairakite.
ERROR: Could not reduce equation to secondary master species, Weeksite.
ERROR: Could not reduce equation to secondary master species, Wollastonite.
ERROR: Could not reduce equation to secondary master species, Xonotlite.
ERROR: Could not reduce equation to secondary master species, Zircon.
ERROR: Could not reduce equation to secondary master species, Zn2SiO4.
ERROR: Could not reduce equation to secondary master species, Zoisite.
ERROR: Calculations terminating due to input errors.
-------------------------------
End of Run after 0.348 Seconds.


Appreciate your help.

[dlparkhurst]

SiO2(aq) is the master species in llnl.dat. All other silica species are written in terms of SiO2, for example SiO2 + 2H2O = H2SiO4-- + 2H+.

When you define H4SiO4 as the master species it invalidates all of the other Si reactions in llnl.dat.

SiO2 and H4SiO4 are essentially equivalent. H4SiO4 is SiO2 plus two waters. Most other databases use H4SiO4(aq) as the master species, but, for these databases, all of the silicate reactions are written with H4SiO4(aq), whereas llnl.dat uses SiO2(aq). I would not try to change the llnl.dat database to use H4SiO4. In addition define H4SiO4 in addition to SiO2, you will have two equivalent species instead of one.

To convert between reactions, you can use the reaction

SiO2 + 2H2O = H4SiO4
log_k 0

So, if you want to use an llnl.dat Si aqueous species or mineral with the wateq4f.dat database, you would use the above equation to remove SiO2 from the llnl.dat reaction and add the new definition in a SOLUTION_SPECIES or PHASES data block. Do the reverse to use a wateq4f.dat reaction with llnl.dat.

[dat]

Thank you for your comment. It was helpful. I have another issue.

What is the difference between the following three codes (Only a part of the full code)?
Because I want to find the optimum answers with the codes. What is the use and meaning of USE EQUILIBRIUM PHASE 1? I am getting different answers from the following codes.

Appreciate your help.

Code: [Select]

SOLUTION 1
    temp      250
    pH        7
    pe        2
    redox     pe
    units     ppm
    density   2.1 calc
    Cl        35500
    Na        23000
    -water    0.25 # kg

GAS_PHASE 1
    -fixed_pressure
    -pressure 110
    -volume 1
    -temperature 40
    CO2(g)    70

INCREMENTAL_REACTIONS True


KINETICS 1
Quartz
    -formula  SiO2  1
    -m        0.0275
    -m0       0.0275
    -parms    4.122
    -tol      1e-08



Code: [Select]

SOLUTION 1
    temp      250
    pH        7
    pe        2
    redox     pe
    units     ppm
    density   2.1 calc
    Cl        35500
    Na        23000
    -water    0.25 # kg

EQUILIBRIUM_PHASES 1
    Anorthite    0 0.00056
    Calcite 0 0.03675
Ankerite 0 0.00395
Quartz 0 0.0275
Albite 0 0.00224
Illite 0 7.5e-06
Clinochlore-14A 0 0.00395


SAVE SOLUTION 1
END


USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1
GAS_PHASE 1
    -fixed_pressure
    -pressure 110
    -volume 1
    -temperature 40
    CO2(g)    70

INCREMENTAL_REACTIONS True


KINETICS 1
Quartz
    -formula  SiO2  1
    -m        0.0275
    -m0       0.0275
    -parms    4.122
    -tol      1e-08



Code: [Select]

SOLUTION 1
    temp      250
    pH        7
    pe        2
    redox     pe
    units     ppm
    density   2.1 calc
    Cl        35500
    Na        23000
    -water    0.25 # kg

EQUILIBRIUM_PHASES 1
    Anorthite    0 0.00056
    Calcite 0 0.03675
Ankerite 0 0.00395
Quartz 0 0.0275
Albite 0 0.00224
Illite 0 7.5e-06
Clinochlore-14A 0 0.00395


SAVE SOLUTION 1
END


USE SOLUTION 1

GAS_PHASE 1
    -fixed_pressure
    -pressure 110
    -volume 1
    -temperature 40
    CO2(g)    70

INCREMENTAL_REACTIONS True


KINETICS 1
Quartz
    -formula  SiO2  1
    -m        0.0275
    -m0       0.0275
    -parms    4.122
    -tol      1e-08



[dlparkhurst]

Calculation 1.

A SOLUTION is defined, which results in an initial solution calculation that distributes species. A second internal calculation determines the number of moles of CO2 in a 1 L volume, given a partial pressure of CO2 of 70 atm and 40 C. Depending on the database, either the ideal gas law is used--PV = nRT, or the Peng-Robinson equation of state is used.
Finally, the second calculation in the output file uses SOLUTION 1, the moles of gas in the gas phase, and runs a kinetic quartz reaction for 1 second (the default time step). The calculation occurs at the temperature and pressure of the solution, 110 atm and 250 C. Equilibrium is established between the gas phase and the solution.

I'm not sure which database and RATES definition you are using. For phreeqc.dat, the Quartz rate (from Rimstidt and, 1980, GCA 44,1683) requires two parameters (the second is a salt correction).

INCREMENTAL_REACTIONS has no effect because there is only one step (1 second) in the reaction calculation.

Calculation 2.

(before the first END) An initial solution calculation is performed for SOLUTION 1. Then solution 1 is equilibrated with the set of minerals in EQUILIBRIUM_PHASES 1, and the reaction solution is saved as solution 1.

Internally, the amount of CO2 in the GAS_PHASE is calculated for 110 atm and 40 C.

(between END and the end of the file) Finally, the saved solution is reacted with EQUILIBRIUM_PHASES 1 (unchanged from initial definition), the gas phase, and with Quartz kinetics reacting for 1 second. This reaction differs from calculation 1 by the inclusion of EQUILIBRIUM_PHASES. You are providing amounts of minerals, which either dissolve to equilibrium or completely dissolve.

INCREMENTAL_REACTIONS has no effect.

Calculation 3.

I think calculation 3 is the same as calculation 2.

[dat]

Thank you for the reply. I am using the llnl database. So it requires only surface area right ? Why we add salt correction?

In the calculation 2 and 3, i am getting different values. The difference is I am adding a USE EQUILIBRIUM PHASE 1. So what could be the reason for that.

Code: [Select]

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.41)*EXP(-90.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-(10^si_qtz))
    190 moles = r * TIME
    200 SAVE moles
     -end
Anorthite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anort = SI("Anorthite")
   30 if (M <= 0 and si_anort < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_anort > 0) then SA = 1e-05  #nucleation
   60 k_acid = (10^-8)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))                 
   70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anort))*ACT("H+")^(1.411)
   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^-12)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_alb))*(ACT("H+")^(0.5))
    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))             
    70 k_neut = 0
    80 k_carb = 0
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-10^(si_calc))*ACT("H+")^(0.5)
    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^(-12.4)*EXP(-46.00e+03/8.314*(1.0/TK-1.0/373.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_ill))*ACT("H+")^(0.1)
    190 moles = r * TIME
    200 SAVE moles
    -end

Clinochlore-14A
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_chlo = SI("Clinochlore-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 = 10^(-12.8)*EXP(-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_chlo))*ACT("H+")^0.5
   190 moles = r * TIME
   200 SAVE moles
-end

Ankerite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_ank = SI("Ankerite")
   30 if (M <= 0 and si_ank < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_ank > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = 10^(-3.5)*EXP(-36.1e+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_ank))*ACT("H+")^0.2
   190 moles = r * TIME
   200 SAVE moles
   -end



SOLUTION 1
    temp      250
    pH        4.41
    pe        2
    redox     pe
    units     ppm
    density   2.1 calc
    Cl        35500
    Na        23000
    -water    0.25 # kg


GAS_PHASE 1
    -fixed_pressure
    -pressure 110
    -volume 1
    -temperature 40
    CO2(g)    70

SAVE SOLUTION 1
SAVE EQUILIBRIUM PHASE 1

END

USE SOLUTION 1
SAVE EQUILIBRIUM PHASE 1

INCREMENTAL_REACTIONS True



KINETICS 1
Quartz
    -formula  SiO2  1
    -m        0.0275
    -m0       0.0275
    -parms    4.122
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.00224
    -m0       0.00224
    -parms    9.5469
    -tol      1e-08
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        0.00056
    -m0       0.00056
    -parms    9.6188
    -tol      1e-08
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        7.5e-06
    -m0       7.5e-06
    -parms    179.5235
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.03675
    -m0       0.03675
    -parms    6.71
    -tol      1e-08
Clinochlore-14A
    -formula  Mg5Al2Si3O10(OH)8  1
    -m        0.00395
    -m0       0.00395
    -parms    678.3
    -tol      1e-08

-steps       604800
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5



[dlparkhurst]

I cannot run your files. I do not have a RATES expression for quartz nor a PHASES definition for Ankerite in my version of llnl.dat.

You need to look carefully at the output file to determine why the results differ.

[dat]

When I run with the same code with USE EQUILIBRIUM PHASE 1, I get the following warnings. I don't get these warnings when I run without  USE EQUILIBRIUM PHASE 1.

And the values are also different. What could be the reason for this? Thank you very much.

Code: [Select]

PHASES

Ankerite
        CaMg0.3Fe0.7(CO3)2 +2.0000 H+  =  + 1.0000 Ca+2 + 0.3Mg+2 +0.7Fe+2 +2.0000 HCO3-
        log_k           2.5135
-delta_H -59.9651 kJ/mol # Calculated enthalpy of reaction Dolomite

Smectite-high-Fe-Mg
        Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12 +8.0000 H+  =  + 0.0250 Ca++ + 0.1000 Na+ + 0.2000 Fe+++ + 0.2000 K+ + 0.5000 Fe++ + 1.1500 Mg++ + 1.2500 Al+++ + 3.5000 SiO2 + 5.0000 H2O
        log_k           17.4200
-delta_H -199.841 kJ/mol # Calculated enthalpy of reaction Smectite-high-Fe-Mg
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.41)*EXP(-90.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-(10^si_qtz))
    190 moles = r * TIME
    200 SAVE moles
     -end
Anorthite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anort = SI("Anorthite")
   30 if (M <= 0 and si_anort < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_anort > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-8)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))                 
   70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anort))*ACT("H+")^(1.411)
   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^(-12)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_alb))*(ACT("H+")^(0.5))
    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))             
    70 k_neut = 0
    80 k_carb = 0
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-10^(si_calc))*ACT("H+")^(0.5)
    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^(-12.4)*EXP(-46.00e+03/8.314*(1.0/TK-1.0/373.15))                 
    70 k_neut = 0
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_ill))*ACT("H+")^(0.1)
    190 moles = r * TIME
    200 SAVE moles
    -end

Clinochlore-14A
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_chlo = SI("Clinochlore-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 = 10^(-12.8)*EXP(-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_chlo))*ACT("H+")^0.5
   190 moles = r * TIME
   200 SAVE moles
-end

Ankerite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_ank = SI("Ankerite")
   30 if (M <= 0 and si_ank < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_ank > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = 10^(-3.5)*EXP(-36.1e+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_ank))*ACT("H+")^0.2
   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 = 10^(-11.31)*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
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))*0
    80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)*0
    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
Magnesite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_magnet = SI("Magnesite")
   30 if (M <= 0 and si_magnet < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_magnet > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-6.38)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.5)                 
   70 k_neut = 10^(-10.78)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_magnet))
   190 moles = r * TIME
   200 SAVE moles
   -end
Siderite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_sid = SI("Siderite")
   30 if (M <= 0 and si_sid < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_sid > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-7.75)*EXP(-48e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.5)                 
   70 k_neut = 10^(-10.78)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_sid))
   190 moles = r * TIME
   200 SAVE moles
   -end
Analcime
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_ana = SI("Analcime")
   30 if (M <= 0 and si_ana < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_ana > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-13.1)*EXP(-58e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.5)                 
   70 k_neut = 10^(-10.78)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_ana))
   190 moles = r * TIME
   200 SAVE moles
   -end
Smectite-high-Fe-Mg
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_sme = SI("Smectite-high-Fe-Mg")
   30 if (M <= 0 and si_sme < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_sme > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-10.98)*EXP(-23.6e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.5)                 
   70 k_neut = 10^(-10.78)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))*0
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_sme))
   190 moles = r * TIME
   200 SAVE moles
   -end


SOLUTION 1
    temp      250
    pH        7
    pe        2
    redox     pe
    units     ppm
    density   2.1 calc
    Cl        35500
    Na        23000
    -water    0.25 # kg

EQUILIBRIUM_PHASES 1
    Anorthite    0 0.00056
    Calcite 0 0.03675
Ankerite 0 0.00395
Quartz 0 0.0275
Albite 0 0.00224
Illite 0 7.5e-06
Clinochlore-14A 0 0.00395


SAVE SOLUTION 1
END


USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1
GAS_PHASE 1
    -fixed_pressure
    -pressure 110
    -volume 1
    -temperature 40
    CO2(g)    70

INCREMENTAL_REACTIONS True


KINETICS 1
Quartz
    -formula  SiO2  1
    -m        0.0275
    -m0       0.0275
    -parms    4.122
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.00224
    -m0       0.00224
    -parms    9.5469
    -tol      1e-08
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        0.00056
    -m0       0.00056
    -parms    9.6188
    -tol      1e-08
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        7.5e-06
    -m0       7.5e-06
    -parms    179.5235
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.03675
    -m0       0.03675
    -parms    6.71
    -tol      1e-08
Clinochlore-14A
    -formula  Mg5Al2Si3O10(OH)8  1
    -m        0.00395
    -m0       0.00395
    -parms    678.3
    -tol      1e-08
Ankerite
    -formula  CaMg0.3Fe0.7(CO3)2  1
    -m        0.00395
    -m0       0.00395
    -parms    6.12
    -tol      1e-08
Analcime
    -formula  Na.96Al.96Si2.04O6:H2O  1
    -m        0
    -m0       0
    -parms    9.24
    -tol      1e-08
K-feldspar
    -formula  KAlSi3O8  1
    -m        0
    -m0       0
    -parms    19.8
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0
    -m0       0
    -parms    298.3
    -tol      1e-08
Magnesite
    -formula  MgCO3  1
    -m        0
    -m0       0
    -parms    2.67
    -tol      1e-08
Siderite
    -formula  FeCO3  1
    -m        0
    -m0       0
    -parms    2.77
    -tol      1e-08
Smectite-high-Fe-Mg
    -formula  Ca.025Na.1K.2Fe.5Fe.2Mg1.15Al1.25Si3.5H2O12 1
    -m        0
    -m0       0
    -parms    1210
    -tol      1e-08
-steps       604800 1209600 3024000
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5


Code: [Select]

Reaction step 1.

WARNING: Element Ca is contained in Ankerite (which has 0.0 mass),
but is not in solution or other phases.
WARNING: Element Ca is contained in Anorthite (which has 0.0 mass),
but is not in solution or other phases.
WARNING: Element Ca is contained in Calcite (which has 0.0 mass),
but is not in solution or other phases.
WARNING: Negative moles in solution 1 for Ca, -1.152683e-02. Recovering...
WARNING: Element Ca is contained in Ankerite (which has 0.0 mass),
but is not in solution or other phases.
WARNING: Element Ca is contained in Anorthite (which has 0.0 mass),
but is not in solution or other phases.
WARNING: Element Ca is contained in Calcite (which has 0.0 mass),
but is not in solution or other phases.
WARNING: Negative moles in solution 1 for Ca, -1.152683e-02. Recovering...
Using solution 1. Solution after simulation 2.
Using pure phase assemblage 1.
Using gas phase 1.
Using kinetics 1.

Kinetics 1.



[dlparkhurst]

The warning is related to the negative moles during the integration. The method goes back and takes smaller steps to fix the problem.

Think about it. In one case you are requiring minerals to dissolve to equilibrium throughout the kinetic reaction, and in the other, only the kinetic reactions occur. The kinetic reactions will have different solution compositions from which the rates are calculated and, of course, the equilibrium phases react when EQUILIBRIUM_PHASES are included, whereas they do not react when not included. Including calcite in EQUILIBRIUM_PHASES and KINETICS does not make sense. The kinetic calcite reaction should not react because calcite is always at equilibrium.

[dat]

Thank you for the clarification.

Actually the warnings appear only when I use the line USE EQUILIBRIUM PHASE 1. The one I didn’t use the line USE EQUILIBRIUM PHASE 1 in kinetics, didn’t show warnings.

That means without using the line USE EQUILIBRIUM PHASE 1 in kinetics, won’t it use the output results of equilibrium phase for kinetics. If not what it will use?

This is hard to understand because of the line USE EQUILIBRIUM PHASE 1.

Thank you.

[dlparkhurst]

In this case, the warning "WARNING: Element Ca is contained in Ankerite (which has 0.0 mass)," occurs when Ankerite was in EQUILIBRIUM_PHASES and (1) there is no Ca in the solution, and (2) no moles of Ankerite in EQUILIBRIUM_PHASES. However, the reason there was no Ca in solution was because of a bad step in the KINETICS integration that removed all of the Ca from solution. Once the integration backed up and took smaller time steps, there was indeed Ca in the solution and no warning message. This situation could occur multiple times during the integration, but as long as there are no "ERROR" messages, the integration eventually completed.

KINETICS evaluates the rate expressions in RATES to determine the rates of reactions and uses a numerical method (either Runge Kutta or CVODE) to integrate those rates over time. No equilibrium phases are required. If you do include a set of equilibrium phases (either with an EQUILIBRIUM_PHASES or USE equilibrium_phases definition), the minerals in equilibrium phases will precipitate or dissolve (provided there are nonzero moles to react) to equilibrium  at every point in the integration (every sub time step). The reaction of the equilibrium minerals to maintain equilibrium will change the water composition, which in turn will affect the rates calculated by the RATES expressions.

The obvious difference in the calculations can be seen in the mole transfers of the equilibrium phases. If you do not include USE equilibrium_phases, there will be no mole transfers. If you include USE equilibrium_phases, then you will have mole transfers of the equilibrium phases. Run just one time step or turn off INCREMENTAL_REACTIONS to make the mole transfers easier to see in the output file.

It does not make sense to put the same mineral in equilibrium phases and in KINETICS, assuming there is some of the equilibrium phase present. The mineral will always be in equilibrium because of the equilibrium phase definition. Assuming the rate expression evaluates to zero when the mineral is in equilibrium (which it should), there will be no mole transfer due to the kinetic reaction. So the kinetic reaction has no effect; all of the reaction occurs because of the equilibrium phase definition.

[dat]

Thank you for the clarification.

Actually, I want the minerals to be equilibrated first and then run a kinetic simulation with CO2.

For this do I have to use the USE EQUILIBRIUM PHASE line? If so what should I put as the mol of minerals in the kinetic phase?

Appreciate your support.

[dlparkhurst]

Adding USE equilibrium_phases during the kinetic reactions will cause the solution to be in equilibrium with those minerals at all points of the kinetic reaction. By not including USE, the saturation indices of those minerals will become nonzero. You must consider whether you expect them to remain at equilibrium or not.

As for moles of the kinetic reactants, those that are present in the formation would have moles present, whereas secondary minerals that are not in the formation would have zero moles. Again, it depends on your conceptual model. The exact moles of a mineral that is present is probably not critical; only a small amount will probably dissolve relative to the moles of mineral in the formation.