Multiple rate blocks conflicting

[DHodkin]

Hello all

I've been working on a script to model kinetic Ca(OH)2 dissolution, CO2 in-gassing and Calcite precipitation from an aqueous solution containing varying initial TIC (fixed with Na2CO3). I have constraints on all these rates, and I have managed to get all of the rate blocks working independent of each other, but when I try to include Ca(OH)2 and Calcite in the same model I get the error message "WARNING: Negative moles in solution for C, -1.#IND00e+00. Recovering..." and the model crashes (I have to quit it with Ctrl+Alt+Del). Are there any reasons why multiple rate blocks would conflict with each other, but work independently?

CO2 in-gassing is defined with a simple linear rate (0.25 mmol/h)
Ca(OH)2 in-gassing is defined by a logarithmic rate (rate=(((11.908*ln(time)+11.237)/100)*M0)/time
The in-built PWP rate was used for Calcite precipitation
Model runs using the Hatches18_no_redox.dat database

RATES
# ----------PHREEQC Calcite rate block-------------------
Calcite
-start
   1 rem   parm(1) = A/V, 1/dm   parm(2) = exponent for m/m0
   10 si_cc = si("Calcite")
   20 if (m <= 0  and si_cc < 0) then goto 200
   30  k1 = 10^(0.198 - 444.0 / (273.16 + tc) )
   40  k2 = 10^(2.84 - 2177.0 / (273.16 + tc) )
   50  if tc <= 25 then k3 = 10^(-5.86 - 317.0 / (273.16 + tc) )
   60  if tc > 25 then k3 = 10^(-1.1 - 1737.0 / (273.16 + tc) )
   70   t = 1
   80   if m0 > 0 then t = m/m0
   90   if t = 0 then t = 1
   100   moles = parm(1) * 0.1 * (t)^parm(2)
   110   moles = moles * (k1 * act("H+") + k2 * act("CO2") + k3 * act("H2O"))
   120   moles = (moles * (1 - 10^(2/3*si_cc)))
   130   moles = moles * time
   140  if (moles > m) then moles = m
   150 if (moles >= 0) then goto 200
   160  temp = tot("Ca")
   170  mc  = tot("C(4)")
   180  if mc < temp then temp = mc
   190  if -moles > temp then moles = -temp
   200 save moles
   210 PRINT moles 
  -end

# ----------Ca(OH)2 rate block-------------------
Ca(OH)2
   -start
    1 si_caoh = si("Ca(OH)2")
      20 if (m <= 0  and si_caoh < 0) then goto 200
    30 rate = (11.908*log(time)+11.237)/100 
    40 moles = rate * M0
    200 SAVE moles
   -end
# ----------Linear CO2 in-gassing--------------------
CO2(g)
   -start
    1 si_co2 = si("CO2(g)")
       20 if (m <= 0  and si_co2 < -3.4) then goto 200
    25 rate = 0.2458/3600000                   #0.25 mmol/hr linear rate
    30 moles = rate * TIME
    200 SAVE moles
   -end
#----------------Starting solution 100ml 0.01M NaOH, 1ppm SrCl2 & 1mmol Na2CO3 --------------------------
solution 1 # pure water
temp 20.0
pH 12.3 charge
water 0.1 kg
units mmol/L
Na 30
C(4) 10
Sr 0.114
Cl 0.228

KINETICS
Ca(OH)2 #Solid Ca(OH)2 added at a concentration of 9mmol per 100ml of solution
      -tol    1e-7   
      -m0 0.009
Calcite #Allow the solution to precipitate Calcite kinetically
      -tol    1e-7
      -m0    0
      -parms   50     0.6
CO2(g) #Allow linear CO2 in-gassing
      -tol    1e-7
         -m0 1
-steps 0.5 1.5 8 50 360 900 1440 2160 2880 3600 4320 minutes
SELECTED_OUTPUT
-file Complete Solution.xls
-reset false
USER_PUNCH
   -headings t Ca Na Sr TIC pH SI_Calcite SI_Aragonite SI_Vaterite SI_Ca(OH)2 SC
10 PUNCH total_time/60/60, TOTMOLE("Ca")*1000, TOTMOLE("Na")*1000, TOTMOLE("Sr")*1000, TOTMOLE("C(4)")*1000, -LA("H+"), SI("calcite"), SI("aragonite"), SI("vaterite"), SI("Ca(OH)2"), SC/1e-3
USER_GRAPH 1
-axis_scale x_axis 0.001 auto auto auto logscale
-axis_scale y_axis auto
-axis_titles "time (h)" "[Ca] mM"
-initial_solutions false
-connect_simulations true
-start
10 GRAPH_X total_time/3600
20 GRAPH_Y TOTMOLE("Ca")*1000
END

[dlparkhurst]

Try using -cvode in KINETICS.

However, I think the more fundamental problem is in your rates. They have strong discontinuities. It would be better if the Ca(OH)2 rate went to zero at equilibrium rather than from a finite rate to zero. Also note that TIME is  an internally kept variable that is a sub-time interval for the integration. TOTAL_TIME is the elapsed time of the simulation. The CO2 ingassing should also go smoothly to zero at pCO2 -3.4 rather than a constant rate discontinuously going to zero at the given pCO2.

[DHodkin]

Thank you for your help, I tried -Cvode in the KINETICS block and the script ran but seemed to skip out the Ca(OH)2 kinetics i.e. [Ca(OH)2] was initially set to ~9 mmol rather than gradually increasing to this over the first 50 minutes.

I think I understand the problem now. The "20 if (m <= 0  and si_caoh < 0) then goto 200" line sets the equilibrium rate to zero but my rate doesn't approach this smoothly causing a discontinuity which is what PHREEQC is having problems computing.

Thank you very much