Reaction-kinetic program failing after 1000 seconds of simulation

[Nitin]

I wrote a code for CO2-Basalt-Water interaction but it is not simulating data even more than one hour. I want to simulate data for around 50 years. So please look into my code. I am using llnl.dat dataset.
Thank you

Code: [Select]

RATES 1
Labradorite
# from Palandri and Kharaka 2004
# experimental condition range T=25-250C, pH=0.5-6

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=5.01E-01
12 E1=42064
13 n1=0.626
20 rem neutral solution parameters
21 a2=1.00E-03
22 E2=45155
30 rem base solution parameters
31 a3=5.62E-10
32 E3=45988
33 n2=-0.782
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("labradorite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles


PHASES
Fix_pe
e- = e-
log_k 0.0
glass
SiTi0.024Al0.358Fe0.188Mg0.281Ca0.264Na0.079K0.008O3.370 + 2.644H+  = 0.358Al+3 + 0.008K+ + 0.079Na+ + 0.264Ca+2 + 0.281Mg+2 + 0.024Ti(OH)4 + 0.171Fe+2 + 0.017Fe+3  + SiO2 + 1.274H2O
           -log_k   -2.60
   -delta_h  30.3 kj/mol
Labradorite
#Na0.4Ca0.6Al1.6Si2.4O8 + 3.2H2O = 0.4Na+ + 0.6Ca+2 + 1.6Al(OH)4- + 2.4SiO2
      Na0.4Ca0.6Al1.6Si2.4O8 = 0.4Na+ + 0.6Ca+2 + 1.6AlO2- + 2.4SiO2
-log_k -20.0702
-analytic -92.306 0.0 -2947.59 33.017 318144 -35.495E-6
END
SOLUTION 1
pH   7.5
Na  1
Ca  0.6
K   0.1
Mg  0.02
Fe  0.0012
Cl  0.3
S   0.1  as  SO4-2
Si   0.2
Al   0.001
Alkalinity   2  as HCO3
use SOLUTION 1
solid_solutions 1 two solid solution
Pigeonite
-comp Enstatite  6.9752
-comp Ferrosilite  5.6352
-comp Wollastonite  5.4496
END

EQUILIBRIUM_PHASES 1
Magnetite 0  0.043
glass  0   7.7
Calcite   0   0
Dolomite  0   0
Magnesite  0   0
Siderite  0   0
CO2(g)   2  10
END
KINETICS 1
Labradorite
 -formula    CaAlSi4O8  1
     -M0 30
     -parms 2000  89
     -steps 15 year in 10 steps
     -cvode true
 INCREMENTAL_REACTIONS true
END
INCREMENTAL_REACTIONS true
USE solution 1
USE equilibrium_phases 1
USE kinetics 1
USER_GRAPH 1
    -headings                rxn  SI("Labradorite") SI("Dolomite")  SI("Calcite")  SI("Siderite") SI("Magnesite") KIN_DELTA("Labradorite")
    -axis_titles            "TIME" "Saturation index"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X total_time
20 GRAPH_Y  SI("Labradorite"), SI("Dolomite"),  SI("Calcite"),  SI("Siderite"),  SI("Magnesite")
30 GRAPH_SY KIN_DELTA("Labradorite")
  -end 
 END
 

[dlparkhurst]

I don't like to admit it, but I think the use of both units and "in 10 steps" does not work properly. I suggest you use seconds for -steps in KINETICS.

Your rate reaches equilibrium in about 10 hours. In that case, you should use EQUILIBRIUM_PHASES for Labradorite. If you want equilibrium to occur on a longer time scale, you need a slower rate.

Code: [Select]

RATES 1
Labradorite
# from Palandri and Kharaka 2004
# experimental condition range T=25-250C, pH=0.5-6

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=5.01E-01
12 E1=42064
13 n1=0.626
20 rem neutral solution parameters
21 a2=1.00E-03
22 E2=45155
30 rem base solution parameters
31 a3=5.62E-10
32 E3=45988
33 n2=-0.782
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("labradorite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles


PHASES
Fix_pe
e- = e-
log_k 0.0
glass
SiTi0.024Al0.358Fe0.188Mg0.281Ca0.264Na0.079K0.008O3.370 + 2.644H+  = 0.358Al+3 + 0.008K+ + 0.079Na+ + 0.264Ca+2 + 0.281Mg+2 + 0.024Ti(OH)4 + 0.171Fe+2 + 0.017Fe+3  + SiO2 + 1.274H2O
           -log_k   -2.60
   -delta_h  30.3 kj/mol
Labradorite
#Na0.4Ca0.6Al1.6Si2.4O8 + 3.2H2O = 0.4Na+ + 0.6Ca+2 + 1.6Al(OH)4- + 2.4SiO2
      Na0.4Ca0.6Al1.6Si2.4O8 = 0.4Na+ + 0.6Ca+2 + 1.6AlO2- + 2.4SiO2
-log_k -20.0702
-analytic -92.306 0.0 -2947.59 33.017 318144 -35.495E-6
END
SOLUTION 1
pH   7.5
Na  1
Ca  0.6
K   0.1
Mg  0.02
Fe  0.0012
Cl  0.3
S   0.1  as  SO4-2
Si   0.2
Al   0.001
Alkalinity   2  as HCO3
use SOLUTION 1
solid_solutions 1 two solid solution
Pigeonite
-comp Enstatite  6.9752
-comp Ferrosilite  5.6352
-comp Wollastonite  5.4496
END

EQUILIBRIUM_PHASES 1
Magnetite 0  0.043
glass  0   7.7
Calcite   0   0
Dolomite  0   0
Magnesite  0   0
Siderite  0   0
CO2(g)   2  10
END
KINETICS 1
Labradorite
 -formula    CaAlSi4O8  1
     -M0 30
     -parms 2000  89
     #-steps 4.725e8 in 10 steps
     -steps 1 3 10 30 100 300 1e3 3e3 1e4 3e4 1e5
     -cvode true
END
INCREMENTAL_REACTIONS true
USE solution 1
USE equilibrium_phases 1
USE kinetics 1
USER_GRAPH 1
    -headings               rxn SI("Labradorite") SI("Dolomite") SI("Calcite") SI("Siderite") \
    SI("Magnesite") Rate("Labradorite")
    -axis_titles            "TIME, hours" "Saturation index" "Rate, mol/s"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X total_time/3600
20 GRAPH_Y  SI("Labradorite"), SI("Dolomite"),  SI("Calcite"),  SI("Siderite"),  SI("Magnesite")
30 GRAPH_SY KIN_DELTA("Labradorite")/KIN_TIME
  -end
    -active                 true
 END
 

[zscowcrofst]

I have couple of questions:
1- Why would you fix pe in the phases. It doesn't have any impact on the results?
2- Can you please cite the references for other parameters i.e. Glass and Labradorite reactions and solid solutions?
Thanks

[dlparkhurst]

Although Fix_pe is defined in the PHASES data block, it is not used in the calculation and could be removed. In general, it is not a good idea to fix the pe. pe should be allowed to adjust to reactions and changes in pH. If fixed, it could cause conditions outside the stability field for water. If you need to control the redox conditions, it is better to fix the partial pressure of O2(g), H2(g), or other redox active gases or phases.

Details on rates, solid solutions, and compositions are simply retained from the original post.

The labradorite rate is stated to come from Palandri and Kharaka 2004 and uses a generalized composition for labradorite.

[zscowcrofst]

Hi David,

I am trying to simulate a situation where the formation water is in equilibrium with host rock including the Glass and Labradorite and then the formation water gets charged with CO2(g) at 100 atm . My understanding is, after charging water with CO2 then it is out of equilibrium and should react kinetically with the minerals in the host rock.

I tried to add glass and labradorite both in first equilibrium, then add CO2(g) in second equilibrium and then let it kinetically react with glass and labradorite but  the kinetics doesn't work because the whole water is in equilibrium with glass and labradorite like this:

Code: [Select]

EQUILIBRIUM_PHASES 1
Magnetite 0  0.043
Glass  0   7.7
Labradorite0 10.0
Calcite   0   0
Dolomite  0   0
Magnesite  0   0
Siderite  0   0
END

EQUILIBRIUM_PHASES 2
CO2(g)   2  10
END

KINETICS 1
Labradorite
    -formula  CaAlSi4O8  1
    -m        30
    -m0       30
    -parms    2000 89
    -tol      1e-08

Glass
    -formula  SiTi0.024Al0.358Fe0.188Mg0.281Ca0.264Na0.079K0.008O3.370  1
    -m        10
    -m0       10
    -parms    500 20
    -tol      1e-08 

I have used these Rates and Phases definitions:

Code: [Select]

RATES 1
Labradorite

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=5.01E-01
12 E1=42064
13 n1=0.626
20 rem neutral solution parameters
21 a2=1.00E-03
22 E2=45155
30 rem base solution parameters
31 a3=5.62E-10
32 E3=45988
33 n2=-0.782
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("labradorite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

Glass
-start
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * 122.566 else S = m0 * ((m/m0)^(2/3)) * 122.566 * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
1000 if (m0<=0) then go to 5000
1001 Aa = 1.08e-4 #mol.m-2.s-1
1003 Ea = 21500   #J.mol-1
1006 R = 8.3144 #J.deg-1.mol-1
1007 ACTI = (ACT ("H+")^3)/(ACT("Al+3"))
1008 n = 1/3
1009 Sig = 1

2000 rplus = Aa * ACTI^n  * exp(-Ea/ (R * Tk)) * S
3000 rate = rplus * (1 - (SR ("Glass")^(1/Sig)))
4000 moles = rate * time
5000 save moles
-end

PHASES
Glass
    SiTi0.024Al0.358Fe0.188Mg0.281Ca0.264Na0.079K0.008O3.370 + 2.644H+ = 0.358Al+3 + 0.264Ca+2 + 0.171Fe+2 + 0.017Fe+3 + 1.274H2O + 0.008K+ + 0.281Mg+2 + 0.079Na+ + SiO2 + 0.024Ti(OH)4
    log_k     -2.6
    delta_h   30.3 kJ
      -analytic   77.82514814711445 0.032450265390183614 -1502.5932036570116 -33.02705435543141 -216815.051931841 -7.454186812457974e-6

Labradorite
    Na0.4Ca0.6Al1.6Si2.4O8 = 1.6AlO2- + 0.6Ca+2 + 0.4Na+ + 2.4SiO2
    log_k     -20.0702
    -analytical_expression -92.306 0 -2947.59 33.017 318144 -3.5495e-05
END

[dlparkhurst]

Please post a complete script (as simple as possible) that demonstrates your problem.

I don't know whether you defined Labradorite in both EQUILIBRIUM_PHASES and KINETICS for your calculation, but you should not put it in both definitions. Use EQUILIBRIUM_PHASES if you want it to be in equilibrium; use KINETICS if you want it to react kinetically.

[zscowcrofst]

Thank you. Here is the input. But at Simulation 5 it iterates over the Glass kinetic calculation and give the following warning:

WARNING: Negative moles in solution 1 for Ti, -2.622257e-13. Recovering...
WARNING: Negative moles in solution 1 for Ti, -1.311127e-13. Recovering...

Code: [Select]

RATES 1
Labradorite
-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=5.01E-01
12 E1=42064
13 n1=0.626
20 rem neutral solution parameters
21 a2=1.00E-03
22 E2=45155
30 rem base solution parameters
31 a3=5.62E-10
32 E3=45988
33 n2=-0.782
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("labradorite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

Glass
-start
2 if (PARM(1) = 0) then  goto 3 else goto 5
3 if PARM(3) = 0 then S = PARM(2) * m * 122.566 else S = m0 * ((m/m0)^(2/3)) * 122.566 * PARM(2)
4 GOTO 1000
5 S = PARM(2)*TOT("water")
1000 if (m0<=0) then go to 5000
1001 Aa = 1.08e-4           #mol.m-2.s-1
1003 Ea = 21500              #J.mol-1
1006 R = 8.3144 #J.deg-1.mol-1
1007 ACTI = (ACT ("H+")^3)/(ACT("Al+3"))
1008 n = 1/3
1009 Sig = 1
2000 rplus = Aa * ACTI^n  * exp(-Ea/ (R * Tk)) * S
3000 rate = rplus * (1 - (SR ("Glass")^(1/Sig)))
4000 moles = rate * time
5000 save moles
-end

PHASES
Glass
    SiTi0.024Al0.358Fe0.188Mg0.281Ca0.264Na0.079K0.008O3.370 + 2.644H+ = 0.358Al+3 + 0.264Ca+2 + 0.171Fe+2 + 0.017Fe+3 + 1.274H2O + 0.008K+ + 0.281Mg+2 + 0.079Na+ + SiO2 + 0.024Ti(OH)4
    log_k     -2.6
    delta_h   30.3 kJ
                -analytic               77.82514814711445 0.032450265390183614 -1502.5932036570116 -33.02705435543141 -216815.051931841 -7.454186812457974e-6

Labradorite
    Na0.4Ca0.6Al1.6Si2.4O8 = 1.6AlO2- + 0.6Ca+2 + 0.4Na+ + 2.4SiO2
    log_k     -20.0702
    -analytical_expression -92.306 0 -2947.59 33.017 318144 -3.5495e-05
END

SOLUTION 1
    temp      25
    pH        7.5
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Al        0.001
    Alkalinity 2 as HCO3
    Ca        0.6
    Cl        0.3
    Fe        0.0012
    K         0.1
    Mg        0.02
    Na        1
    S         0.1 as SO4-2
    Si        0.2
    -water    1 # kg

use SOLUTION 1

SOLID_SOLUTIONS 1 two solid solution
    Pigeonite
        -comp Enstatite 6.9752
        -comp Ferrosilite 5.6352
        -comp Wollastonite 5.4496
    Augite
        -comp Enstatite 4.5
        -comp Ferrosilite 5.6
        -comp Wollastonite 10.5
END

EQUILIBRIUM_PHASES 1
    Calcite   0 0
    Dolomite  0 0
    Glass     0 10
    Labradorite 0 10
    Magnesite 0 0
    Magnetite 0 0.043
    Siderite  0 0
END

EQUILIBRIUM_PHASES 2
    CO2(g)    2.5 10
END

USE solution 1

KINETICS 1
Labradorite
    -formula  CaAlSi4O8  1
    -m        30
    -m0       30
    -parms    2000 89
    -tol      1e-08

Glass
    -formula  SiTi0.024Al0.358Fe0.188Mg0.281Ca0.264Na0.079K0.008O3.370  1
    -m        20
    -m0       20
    -parms    1000 50
    -tol      1e-08
-steps       1 3 10 30 100 300 1000 3000 10000 30000 100000
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
END

INCREMENTAL_REACTIONS True

USE equilibrium_phases 2

USE kinetics 1

USE solution 1

USER_GRAPH 1
    -headings               rxn ("Glass") ("Labradorite") ("Dolomite") ("Calcite") ("Siderite") ("Magnesite")
    -axis_titles            "TIME, hours" "Saturation index" ""
    -chart_title            "Saturation Index"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start

10 GRAPH_X total_time/3600
20 GRAPH_Y  SI("Glass"),SI("Labradorite"), SI("Dolomite"),  SI("Calcite"),  SI("Siderite"),  SI("Magnesite")
  -end
    -active                 true
 
USER_GRAPH 2
    -headings               rxn ("Glass") ("Labradorite") ("Dolomite") ("Calcite") ("Siderite") ("Magnesite")
    -axis_titles            "Time [d]" "Rate [mol/s]" ""
    -chart_title            "Rates"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start

10 GRAPH_X total_time/86400
20 GRAPH_Y  KIN_DELTA("Glass"),KIN_DELTA("Labradorite"), KIN_DELTA("Dolomite"),  KIN_DELTA("Calcite"),  KIN_DELTA("Siderite"),  KIN_DELTA("Magnesite")
  -end
    -active                 true

USER_GRAPH 3
    -headings               rxn ("Ca") ("Si") ("Al") ("C4") ("Fe(2)") ("Ti") ("pH")
    -axis_titles            "Time [Day]" "Concentration [ppm]" "pH [-]"
    -chart_title            "Concentration & pH"
    -axis_scale x_axis      auto auto auto auto log
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start

10 GRAPH_X total_time/86400
30 GRAPH_Y TOT("Ca"), TOT("Si"), TOT("Al"), TOT("C4"),TOT("Fe(2)"),TOT("Ti")
40 GRAPH_SY -LA("H+")
  -end
    -active                 true
END

[dlparkhurst]

That is a WARNING, not an ERROR. As long as you don't get an error, the integration method recovers by running the time interval with smaller sub time steps.

The results should be correct.

[zscowcrofst]

Hi,
Thank you for your feedback.

If I remove the Glass from KINETIC it runs pretty fast in few second but once I add the Glass it takes for very long time. What is could be the reason? Is the reaction rate too fast? Or stoichiometry is not right? Is there anything in ODE Method and Steps that can speed up the simulation?

Thanks

[dlparkhurst]

I think the trace concentrations of Ti and Fe in the glass cause problems.

Rather than kinetics, I think you should use EQUILIBRIUM_PHASES to investigate your system. If you run to equilibrium, you will find that you are way supersaturated with many secondary minerals. You will also see extremely small concentrations of Fe or Ti depending on secondary minerals that you include. So, use equilibrium phases consider the end point that you want your kinetic reaction to tend toward. Consider whether you want glass to dissolve (negative tranfer), or form (positive transfer).


I think your original problem, and the problem that will recur if you go back to kinetic calculations, is that the integration grinds to a halt when concentrations become small ( < 1e-12?). The integrator is trying to remove reactants, but it must take micro time steps to avoid negative concentrations.