Infinite iteration and something weird happening with TIME in RATES block

[AnotherSLGuy]

Hello,

I am new to using PHREEQC and I am trying to calculate the gypsum precipitation rate using RATES and KINETICS.

Accordingly, I wrote the following code for gypsum precipitation based on literature. It has some placeholder values for some of the variables for the time being, but the algebraic expressions are in their final form.

When executed, the rate block goes into an infinite loop. While trying to find the issue, I tried printing intermediate outputs to see whether there are any abnormalities, and it seems TIME is acting weird and I have not clue as to why. (I have included a sample output at the end.)

I would greatly appreciate if I can get some pointers as to what is happening and how to avoid it in the future.

Thank you.

Code: [Select]

SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Ca        0.1
    S(6)      0.1
    -water    1 # kg


RATES
    Gypsum
-start
 10 si_gyp = SI("Gypsum")
 20 IF(M <= 0 AND si_gyp < 0) THEN GOTO 250
 30 b_het = 8.05*10^-19
 40 log_Ksp_Gyp =(-1.6201991*10^3)-(2.5723285*10^(-1)*TK)+(8.9149708*10^(4)/TK)+(5.8737969*10^2*LOG10(TK))-(5.3472954*10^6/(TK^2))
 50 C_s = ((10^log_Ksp_Gyp)/(GAMMA("Ca+2")+GAMMA("SO4-2")))^0.5
 60 beta = 16.755
 70 V_m = 7.469*10^-5
 80 zigma = 12.2*10^-3
 90 NA = 6.022*10^23
100 f_teta = 5.17*10^-1
110 R = 8.314472
120 Omega = SR("Calcite")
130 k1 = 1.04*10^-8
140 mt = M0*172.17
150 Nt = M0*6.022*10^23
160 p = 2/3
170 m_sol = MOL("H2O")*18.01528*10^-3
180 f_DeltaG = ((Omega^0.5)-1)^2
190 M_nuc = -2*beta*(V_m/NA)^2*(zigma)^3/(f_DeltaG/NA)^-3
200 Js_Tot = b_het*C_s*exp(-(beta*V_m^2*zigma^3*NA*f_teta)/((R*TK)^3*LOG(Omega)^2))
210 Rate_het = k1*(mt/Nt)^p*Nt*f_DeltaG
220 Rate = (Js_Tot*M_nuc)+(Rate_het/m_sol)
230 moles = (Rate*TIME)+ moles
240 Print TIME
250 SAVE moles
-end

KINETICS 1
Gypsum
    -formula  CaSO4:2H2O  1
    -m        0.005
    -m0       0.001 
    -tol      1e-08
-steps       864000 in 20 steps # seconds
-step_divide 100
-runge_kutta 6
-bad_step_max 500



USER_GRAPH 1 Molarity variation
    -headings               _time_ Gypsum_KIN
    -axis_titles            "Time in days" "Moles/kg of water" "SI"
    -chart_title            "Gypsum variation"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600/24
20 GRAPH_Y KIN("Gypsum")
  -end
    -active                 true
END

SAMPLE OUTPUT

-----------------------------------------
Beginning of batch-reaction calculations.
-----------------------------------------

Reaction step 1.

        4320
  2.5732e+00
  2.5732e+00
  2.5732e+00
  2.5732e+00
  2.5732e+00
  2.5732e+00
  1.0293e+01
  1.8972e+00
  1.8972e+00
  1.8972e+00
  1.8972e+00
  1.8972e+00
  1.8972e+00
  7.5890e+00
  1.7345e+00
  1.7345e+00
  1.7345e+00
  1.7345e+00
  1.7345e+00
  1.7345e+00
  6.9380e+00
  1.6830e+00
  1.6830e+00
  1.6830e+00
  1.6830e+00
  1.6830e+00
  1.6830e+00
  6.7318e+00
  1.6653e+00
  1.6653e+00
  1.6653e+00
  1.6653e+00
  1.6653e+00
  1.6653e+00
  6.6611e+00
  1.6591e+00
  1.6591e+00
  1.6591e+00
  1.6591e+00
  1.6591e+00
  1.6591e+00
  6.6362e+00
  1.6568e+00
  1.6568e+00
  1.6568e+00
  1.6568e+00
  1.6568e+00
  1.6568e+00
  6.6274e+00
  1.6561e+00
  1.6561e+00
  1.6561e+00
  1.6561e+00
  1.6561e+00
  1.6561e+00
  6.6242e+00
  1.6558e+00
  1.6558e+00
  1.6558e+00
  1.6558e+00
  1.6558e+00
  1.6558e+00
  6.6231e+00
  1.6557e+00
  1.6557e+00
  1.6557e+00
  1.6557e+00
  1.6557e+00
  1.6557e+00
  6.6227e+00
  1.6556e+00
  1.6556e+00
  1.6556e+00

[dlparkhurst]

If the formula has positive coefficients (-formula CaSO4 1 H2O 2), then when the saturation ratio is greater than 1, the value for SAVE should be negative. Conversely, when SR < 1, SAVE should be positive. Your rate does not follow this requirement, and continues to add gypsum even when supersaturated.

A factor like (1 - SR("Gypsum")) is often used to provide the appropriate sign for the reaction. I'm not going to try to decipher your use of Omega.

I also think your rate is very fast. Try running for 1 second and use an arbitrary small rate constant in SAVE to slow down the reaction until you determine appropriate parameters.

[AnotherSLGuy]

Thank you very much for the response.

Regarding the sign of the rate, I manually flipped the sign of "Rate" and it starts behaving somewhat acceptably. Is this due to how the RATES data block is integrated using Runge-Kutta method? (I currently do not have any knowledge of that method TBH.)

Omega is a bit ambiguous at the moment. The literature I used defines it as "the degree of supersaturation of the solution" and does not provide any expression. I suppose it is a common entity in chemistry, but as my background is civil engineering, I did not know what it is. Based on what I read so far on the internet, SR is equivalent of the degree of supersaturation. But there is a good chance for it to be wrong.

I just tried manually slowing down and it does work. 

[dlparkhurst]

The sign convention is relative to the solution. If the product of the quantity SAVEd times a stoichiometric coefficient in -formula is positive, then the reactant is added to the solution. If the quantity SAVEd times a stoichiometric coefficient is negative, then the reactant is removed from solution.

These two add CaSO4 to solution:

Code: [Select]

SAVE 1
-formula CaSO4 1

and

SAVE -1
-formula CaSO4 -1

Omega is normally the same as the saturation ratio (the Basic function SR).