Kinetics of minerals mixture

[Agabamud]

I'm trying to get to grips with the KINETICS and RATES (using sit.dat) and tried to perform a simple dissolution with various of minerals ( Calcite, Portlandite, Hydroxyapatite & Magnesite) in a sequence.

First I performed each dissolution individually. There are no issues with the Calcite and Magnesite dissolutions, though when I repeat the same for Hydroxyapatite and Portlandite, I get the same type of error:

WARNING: Negative moles in solution 1 for Ca, -2.352461e-13. Recovering...

and the system will continue trying to simulate endlessly.

Then I tried to perform the same simulation but with the mineral dissolutions one after the other in steps. In this case I get a different error but of a similar type:

WARNING: Negative moles in solution 3 for P, -1.088065e-10. Recovering...

and once again the system will continue trying to simulate endlessly.

Below is the script I am using:

Code: [Select]

SELECTED_OUTPUT 1
    -file                 PHREEQC_VT_KINETICS_output.sel
    -simulation           true
    -state                true
    -solution             true
    -time                 true
    -step                 true
    -pH                   true
    -reaction             true
    -temperature          true
    -water                true
    -totals               Ca  Mg  C  P
    -molalities           CO3-2  Ca+2  Mg+2  OH-
                          PO4-3
    -equilibrium_phases   Calcite  Hydroxyapatite  Magnesite  Portlandite
    -saturation_indices   CO2(g)  H2O(g)  Portlandite  Mg3(PO4)2:8H2O(s)
                          Mg3(PO4)2(cr)  Mg(HPO4):3H2O(s)  Magnesite(nat)  Magnesite(syn)
                          Magnesite
    -kinetic_reactants    Calcite  Hydroxyapatite  Magnesite  Portlandite
PHASES
Calcite
    Ca1C1O3 + H+ = Ca+2 + HCO3-
    -analytical_expression 4764.488 1.443789 -196549.5 -1855.352 8937845 -0.0004900053
    -Vm       36.89 cm3/mol
Portlandite
    Ca(OH)2 + 2H+ = Ca+2 + 2H2O
    log_k     -22.81
    delta_h   130.078 kJ
Hydroxyapatite
    Ca5P3O13H1 + 4H+ = 5Ca+2 + H2O + 3HPO4-2
    -analytical_expression 24302.76 7.482048 -989661.7 -9493.973 44065860 -0.002584642
    -Vm       159.6 cm3/mol
Magnesite
    Mg1C1O3 + H+ = HCO3- + Mg+2
    -analytical_expression 2820.078 0.9099383 -110094.2 -1109.264 4816883 -0.0003275349
    -Vm       28.03 cm3/mol
RATES
    Calcite
-start
 11 a1=0
 12 E1=0
 13 n1=0
 21 a2=6.59E+04
 22 E2=66000
 31 a3=1.04E+09
 32 E3=67000
 33 n2=1.6
 40 SR_mineral=SR("calcite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("HCO3-")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
    Magnesite
-start
 11 a1=1.39E-4
 12 E1=14400
 13 n1=1.000
 21 a2=5.99E-6
 22 E2=23500
 31 a3=6.03E+05
 32 E3=62800
 33 n2=1.000
 40 SR_mineral=SR("magnesite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*SR("CO2(g)")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral^4)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
    Portlandite
-start
 11 a1=1.10E+10
 12 E1=75000
 13 n1=0.600
 21 a2=3.04E+05
 22 E2=75000
 31 a3=0
 32 E3=0
 33 n2=0
 40 SR_mineral=SR("Portlandite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
    Hydroxyapatite
-start
 11 a1=5.12E-05
 12 E1=0
 13 n1=0.171
 21 a2=1E-06
 22 E2=0
 31 a3=0
 32 E3=0
 33 n2=0
 40 SR_mineral=SR("Hydroxyapatite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
END
SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    0.03 # kg
END
USE solution 1
KINETICS 1
Calcite
    -formula  CaCO3  1
    -m        0.015
    -m0       0
    -parms    3 1
    -tol      1e-08
-steps       120 in 4 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
INCREMENTAL_REACTIONS True
SAVE solution 2
END
USE solution 2
KINETICS 2
Magnesite
    -formula  Mg1C1O3  1
    -m        0.018
    -m0       0
    -parms    25 1
    -tol      1e-08
-steps       120 in 4 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
INCREMENTAL_REACTIONS True
SAVE solution 3
END
USE solution 3
KINETICS 2
Hydroxyapatite
    -formula  Ca5P3O13H1  1
    -m        0.003
    -m0       0
    -parms    126 1
    -tol      1e-08
-steps       120 in 4 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
INCREMENTAL_REACTIONS True
SAVE solution 4
END
USE solution 4
KINETICS 2
Portlandite
    -formula  Ca(OH)2  1
    -m        0.02
    -m0       0
    -parms    19 1
    -tol      1e-08
-steps       120 in 4 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
INCREMENTAL_REACTIONS True
SAVE solution 5
END




Could anyone offer any support? I don't understand why it works in certain situation and not in others.

Thanks!

[dlparkhurst]

I think there are a few issues.

First, I think the log K and delta H for portlandite have the wrong signs.

Normally, M0 is equal to M in the KINETICS definitions. As M decreases, the M/M0 factor is intended to account for changes in surface area. If M0 is not defined, it is set equal to M.

Your definition of Hydroxyapatite has a very low solubility, which causes very small concentrations of P. Numerically, it becomes difficult to integrate the rate when the concentration of P is 1e-10 or less. I see three options, (1) you could move hydroxyapatite to equilibrium phases, which is much easier to solve. With concentrations on the order of 1e-10, the mole transfer of hydroxyapatite is negligible, so the exact rate is not important, (2) decrease the tolerance for Hydroxyapatite in KINETICS to 1e-12. The integration may still be slow, but I think it will have a better chance of succeeding. You may also want to play with the convergence tolerance with KNOBS if you make the integration tolerance small. (3) Use a more soluble definition for Hydroxyapatite. The range of equilibrium constants in different databases is large, so you have a choice for the stability of the mineral. I would probably use option 1, but would consider the appropriate log K.

Use the implicit method -cvode in KINETICS. It will do better when there is a large difference in rates among the minerals.

Use 1 kg water. PHREEQC works best with approximately 1 kg water in the solution. If necessary, adjust the moles of kinetic reactants accordingly.

Code: [Select]

PHASES
Calcite
    Ca1C1O3 + H+ = Ca+2 + HCO3-
    -analytical_expression 4764.488 1.443789 -196549.5 -1855.352 8937845 -0.0004900053
    -Vm       36.89 cm3/mol
Portlandite
    Ca(OH)2 + 2H+ = Ca+2 + 2H2O
#    log_k     -22.81
#    delta_h   130.078 kJ
    log_k     22.81
    delta_h   -130.078 kJ
Hydroxyapatite
    Ca5P3O13H1 + 4H+ = 5Ca+2 + H2O + 3HPO4-2
    -analytical_expression 24302.76 7.482048 -989661.7 -9493.973 44065860 -0.002584642
    -Vm       159.6 cm3/mol
Magnesite
    Mg1C1O3 + H+ = HCO3- + Mg+2
    -analytical_expression 2820.078 0.9099383 -110094.2 -1109.264 4816883 -0.0003275349
    -Vm       28.03 cm3/mol
RATES
    Calcite
-start
 11 a1=0
 12 E1=0
 13 n1=0
 21 a2=6.59E+04
 22 E2=66000
 31 a3=1.04E+09
 32 E3=67000
 33 n2=1.6
 40 SR_mineral=SR("calcite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("HCO3-")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
#210 PRINT "Calcite: ", moles, M, rate, TIME, SA
-end
    Magnesite
-start
 11 a1=1.39E-4
 12 E1=14400
 13 n1=1.000
 21 a2=5.99E-6
 22 E2=23500
 31 a3=6.03E+05
 32 E3=62800
 33 n2=1.000
 40 SR_mineral=SR("magnesite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*SR("CO2(g)")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral^4)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
#210 PRINT "Magnesite: ", moles, M, rate, TIME, SA
-end
    Portlandite
-start
 11 a1=1.10E+10
 12 E1=75000
 13 n1=0.600
 21 a2=3.04E+05
 22 E2=75000
 31 a3=0
 32 E3=0
 33 n2=0
 40 SR_mineral=SR("Portlandite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
#210 PRINT "Portlandite: ", moles, M, rate, TIME, SA
-end
    Hydroxyapatite
-start
 11 a1=5.12E-05
 12 E1=0
 13 n1=0.171
 21 a2=1E-06
 22 E2=0
 31 a3=0
 32 E3=0
 33 n2=0
 40 SR_mineral=SR("Hydroxyapatite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
#210 PRINT "Hydroxyapatite: ", moles , M, rate, TIME, SA
-end
END
SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
   # -water    0.03 # kg
END
#KNOBS
#-conv 1e-12
USE solution 1
#EQUILIBRIUM_PHASES
#Hydroxyapatite 0 0.003
KINETICS 1
Calcite
    -formula  CaCO3  1
    -m        0.015
#    -m0       0
    -parms    3 1
    -tol      1e-08
Magnesite
    -formula  Mg1C1O3  1
    -m        0.018
#    -m0       0
    -parms    25 1
    -tol      1e-08
Hydroxyapatite
    -formula  Ca5P3O13H1  1
    -m        0.003
#    -m0       0.003
    -parms    126 1
    -tol      1e-12
Portlandite
    -formula  Ca(OH)2  1
    -m        0.02
#    -m0       0
    -parms    19 1
    -tol      1e-08
-cvode
-step 1e4 in 10 #120
INCREMENTAL_REACTIONS
END



[Agabamud]

Thanks for your support, I have made the changes you suggested and the script is working well.

[Agabamud]

@dlparkhurst I have a few follow up questions. As mentioned earlier the script with the corrections you suggested worked well. First to have a better understanding of the KINETIC feature. Is it that with each time step there is the defined amount of moles (m) that is added to the solution?

Then I changed the order of the kinetic tasks, so literally swapping the Hydroxyapatite and Portlandite tasks, as I would like to see if addition sequence affects the final mixture composition and parameters (only changing the time steps to have the last one 12 hrs and all the rest to 120s) The script is as seen below:

Code: [Select]

SELECTED_OUTPUT 1
    -file                 PHREEQC_VT_KINETICS_output(15).sel
    -simulation           true
    -state                true
    -solution             true
    -time                 true
    -step                 true
    -pH                   true
    -reaction             true
    -temperature          true
    -water                true
    -totals               Ca  Mg  C  P
    -molalities           CO3-2  Ca+2  Mg+2  OH-
                          PO4-3
    -equilibrium_phases   Calcite  Hydroxyapatite  Magnesite  Portlandite
    -saturation_indices   Aragonite  Brucite  CaCO3:H2O(s)  Calcite
                          H2O(g)  Hydroxyapatite  Portlandite  Vaterite
                          Artinite  Brushite  Ca(HPO4)(s)  Ca3(PO4)2(alfa)
                          Ca4H(PO4)3:2.5H2O(s)  CaCO3:H2O(s)  Calcite  CaMg3(CO3)4(s)
                          Dolomite  Magnesite  Periclase  Mg5(CO3)4(OH)2:4H2O(s)
                          Nesquehonite  Mg3(PO4)2:8H2O(s)  Magnesite(nat)  Magnesite(syn)
    -kinetic_reactants    Calcite  Hydroxyapatite  Magnesite  Portlandite
PHASES
Calcite
    Ca1C1O3 + H+ = Ca+2 + HCO3-
    -analytical_expression 4764.488 1.443789 -196549.5 -1855.352 8937845 -0.0004900053
    -Vm       36.89 cm3/mol
Portlandite
    Ca(OH)2 + 2H+ = Ca+2 + 2H2O
    log_k     22.81
    delta_h   -130.078 kJ
    -analytical_expression 0.0213461 0 6794.44 0 0 0
Hydroxyapatite
    Ca5(OH)(PO4)3 + 7H+ = 5Ca+2 + 3H2(PO4)- + H2O
    log_k     14.35
    delta_h   -178.487 kJ
    -analytical_expression -16.9195 0 9323.01 0 0 0
Magnesite
    Mg1C1O3 + H+ = HCO3- + Mg+2
    -analytical_expression 2820.078 0.9099383 -110094.2 -1109.264 4816883 -0.0003275349
    -Vm       28.03 cm3/mol
RATES
    Calcite
-start
 11 a1=0
 12 E1=0
 13 n1=0
 21 a2=6.59E+04
 22 E2=66000
 31 a3=1.04E+09
 32 E3=67000
 33 n2=1.6
 40 SR_mineral=SR("calcite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("HCO3-")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
    Magnesite
-start
 11 a1=1.39E-4
 12 E1=14400
 13 n1=1.000
 21 a2=5.99E-6
 22 E2=23500
 31 a3=6.03E+05
 32 E3=62800
 33 n2=1.000
 40 SR_mineral=SR("magnesite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*SR("CO2(g)")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral^4)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
    Portlandite
-start
 11 a1=1.10E+10
 12 E1=75000
 13 n1=0.600
 21 a2=3.04E+05
 22 E2=75000
 31 a3=0
 32 E3=0
 33 n2=0
 40 SR_mineral=SR("Portlandite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
    Hydroxyapatite
-start
 11 a1=5.12E-05
 12 E1=0
 13 n1=0.171
 21 a2=1E-06
 22 E2=0
 31 a3=0
 32 E3=0
 33 n2=0
 40 SR_mineral=SR("Hydroxyapatite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
115 if (moles>M) then moles=M
200 save moles
-end
USER_GRAPH 2 Rate graph
    -axis_titles            "Total time (seconds)" "Rate" ""
    -chart_title            "Evolution of rates"
    -axis_scale x_axis      auto 120 auto auto
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y KIN_DELTA("Calcite")/KIN_TIME
  -end
    -active                 true
USER_GRAPH 1 pH trend
    -axis_titles            "time (hours)" "pH" ""
    -chart_title            "pH evolution"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y -LA("H+")
  -end
    -active                 true
END
SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
END
USE solution 1
KINETICS 1 Calcite addition
Calcite
    -formula  CaCO3  1
    -m        1.499
    -m0       1.499
    -parms    3 1
    -tol      1e-008
-steps       120 in 4 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
INCREMENTAL_REACTIONS True
SAVE solution 2
END
USE solution 2
KINETICS 2 Magnesite addition
Magnesite
    -formula  Mg1C1O3  1
    -m        1.779
    -m0       1.779
    -parms    25 1
    -tol      1e-008
-steps       120 in 4 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
INCREMENTAL_REACTIONS True
SAVE solution 3
END
USE solution 3
KINETICS 4 Portlandite addition
Portlandite
    -formula  Ca(OH)2  1
    -m        2.024
    -m0       2.024
    -parms    19 1
    -tol      1e-08
-steps       120 in 4 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
INCREMENTAL_REACTIONS True
SAVE solution 4
END
USE solution 4
KINETICS 3 Hydroxyapatite addition
Hydroxyapatite
    -formula  Ca5P3O13H1  1
    -m        0.299
    -m0       0.299
    -parms    126 1
    -tol      1e-012
-steps       43200 in 12 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
INCREMENTAL_REACTIONS True
SAVE solution 5
END




Though when I run this script I run into a problem that we had at the start of this where the system says there are negative moles for P. In the end the system runs the simulation but I've noticed that the results are constant for the for all the time steps in the last simulation (which is after the error).

This is in contrast to the simulation I ran before I swapped around the hydroxyapatite and Portlandite KINETICS data block where the various parameters were evolving with every time step. I wanted to know if the constant results are a symptom of the simulation not functioning correctly?

[dlparkhurst]

No. The amount of reactant added is not M; the amount is calculated by the program defined in RATES. During the integration, the amounts defined by the SAVE statements are added during sub time steps. The SAVE amount is the amount of reaction that occurs over a sub time step determined internally by the variable TIME. The Runge-Kutta integration uses multiple sub time steps to estimate the total amount of reaction that occurs over the time step defined in KINETICS.

The WARNING about negative moles is not a critical error. It simply means that the sub time steps were too large, resulting in a poor integration. The method automatically resets and integrates with a smaller time step to avoid the negative concentrations. If the method fails, there will be an ERROR message and the calculation will not be completed.

I'm not sure about the constant results, but logically, the sequential calculation is different from the calculation where all of the reactions are occurring simultaneously. In the simultaneous calculation, it should not matter the order that the reactions are defined in KINETICS. In the sequential kinetics calculations, the order probably does matter. The RATES definitions are dependent on solution composition, so when you switch the order of reactions, you are changing the solution compositions for the kinetic reactions. I expect that the different solution compositions will result in different rates and different solution compositions at the end of each KINETICS calculation.


Unless there is an ERROR message, I think the integrations are functioning correctly, and any differences are due to differences in the sequence of calculations, which represent different conceptual models.

[Agabamud]

@dlparkhurst Thanks for the information, I really appreciate it.

Follow up questions for clarity. Therefore to know how much was dosed each time step one must refer to the Total Moles and the Delta Moles in the output file?

I'm guessing the only way to define how much is actually reacting is to use EQUILIBRIUM PHASES and input the desired moles there?

You were right about the integrations functioning correctly even when there are errors. It appears that the sequencing affects results in rates, molarities, saturation indices etc. etc. sometimes producing minimal or no change over time.

Just out of curiosity what function would you use on PHREEQC to most accurately simulate a simple water and mineral(s) interaction. The idea being that you add the minerals one after the other, allow for mixing, monitor the effect of each mineral addition and then hopefully being able to monitor the final mixture evolution over time? I have tried using MIX, REACTION and now the KINETICS & RATES data blocks, and with each one there are limitations. Seeing as you have more experience with the software I would like to know what you would do.

[dlparkhurst]

Regarding kinetic simulations, the total moles and delta moles are printed in the output file. If you use INCREMENTAL_REACTIONS true and perform multiple reaction steps, the delta moles is the change from the last time step; if you use INCREMENTAL_REACTIONS false the delta moles is from the initial time of the simulation.

You can retrieve the values in USER_PRINT and USER_PUNCH scripts with the functions KIN and KIN_DELTA. The function KIN will return the current moles of a kinetic reactant; the function KIN_DELTA will return the change in moles (possibly subject to the choice of INCREMENTAL_REACTIONS described above.

The total amount of a KINETICS reactant is given by the identifier -M. The value of M (also equal to KIN()) will vary with time. When a SAVEd value (from the RATES definition) is positive, M decreases; when a SAVEd value is negative, M increases. The value of M may decrease until it reaches zero, at which point the reaction will stop (unless negative values for SAVE are encountered).

With kinetics, I think reactions occurring simultaneously is more common than a sequence of reactions. If the reactions are occurring simultaneously, I would not model them as a sequence because I think the results would vary depending on the sequence you used. I can't say whether REACTION, KINETICS, or EQUILIBRIUM_PHASES suits your purposes. I think you understand how each works, so you will have to decide which best represents your conceptual model. There are more assumptions involved in KINETICS because of the rate expression than in REACTION or EQUILIBRIUM_PHASES.

[Agabamud]

@dlparkhurst Hi, I've noticed that the KINETICS data block unlike the other data blocks on phreeqc (REACTION, EQUILIBRIUM PHASES) doesn't seem to be highly reactive to changes in concentration. They are more responsive to changes in temperature. Is this mineral dependent or just a feature of the model coming from the conditions of actual geological reactions? or even further than that could it purely depend on the rate scripts that have been defined?

Then I have a 2nd question related to how PHREEQC fundamentally works to generate information. Ill take the example of a reaction involving a mineral that are not defined as phases. So is carried out using the REACTION data block. Am I right in thinking that the model breaks down the reacting species into ions and then will produce data for each (possible) combination of ions? If yes, then how does the model choose between different combinations(for example H+ with either Cl- or O-2) ?

I ask this because I looked through the database and there isn't much information provided for the ionic species in the SOLUTION SPECIES and SOLUTION MASTER SPECIES outside of the parameters used for temperature dependence (analytical equation & delta h) and log K.
I would have expected that enthalpies would be used for this but the enthalpy of formation is not used by the model and the delta_h is used "to determine the temperature dependence of the equilibrium constant".
What am I missing to understand how PHREEQC works?

[dlparkhurst]

The amount of reaction using KINETICS depends on the RATES definition and the solution composition. The rate of dissolution for a mineral far from equilibrium will depend on the rate expression, which could be on the order of hours or millennia. Rates would slow as a mineral nears equilibrium.

Excluding KINETICS, PHREEQC solves a set of nonlinear equations. For a SOLUTION plus REACTION, the set of equations includes mole-balance equations for H, O, and charge, plus one mole-balance equation for every other element in solution. In addition, there are mass-action equations for each aqueous species, equations like K = [HCO3-]/([H+][CO3-2]). There are additional equations for activity coefficients, ionic strength, activity of water and others. A complete description of the equations and numerical method is given in the manual for Version 2 (1999), which is distributed in the doc directory of the installation directory.

SOLUTION_SPECIES defines the reaction for each aqueous species, from which the mass-action expressions are derived. There must also be a definition of the log K for the reaction, which could be just the 25 C log K, Log K plus enthalpy of reaction to calculate temperature dependence, or an analytical expression to calculate log K as a function of temperature. In addition there may be parameters for activity coefficients, diffusion coefficient, and molar volume. These are all the data necessary to formulate the equations for a given chemical system.