The rate equation with "nucleation" does not work at m0 = 0

[MichaelZ20]

David, could you please help?
In the attached input file, the rate equation with "nucleation" does not work only at m0 = 0.

Code: [Select]

SOLUTION 1
    pH         7   charge
    density   1   calculate
    Cl         1
    Na         1
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Calcite   0   10
    CO2(g)     -1   
SAVE SOLUTION 1
END

USE SOLUTION 1

EQUILIBRIUM_PHASES
CO2(g) -3.4

INCREMENTAL_REACTIONS True

RATES 
Calcite
-start
1 dpmet = ((1 / TK) - (1 / 298.15)) / 8.3144
2 SRmin = SR("Calcite")
10 moles = 0
20 If (m <= 0) and (SRmin < 1) Then GoTo 220
30 S = 0.7     # average BET; suggested value in m2/g
40 Mm = 100.1  # molar mass in g/mol
50 If (SRmin > 1) Then GoTo 130
130 knu = 0.000000148 * exp((-63500) * dpmet)
140 k1 = 0.00658 * exp((-56100) * dpmet) * (ACT("HCO3-") ^ 1.982)
150 kpre = (-1) * (knu + k1)
160 theta = 0.5
170 eta = 2
#180 If (m <= 0) then GoTo 220
185 if (m = 0 and SRmin > 1) then rate = 1e-9           # NUCLEATION
#186 GoTo 190
187 if (m > 0 and SRmin > 1) then GOTO 190
190 rate = S * m * Mm * kpre * (ABS(1 - SRmin ^ theta) ^ eta)
#200 GoTo 210
210 moles = rate * Time
220 Save moles
-end   
   
KINETICS 1
Calcite
    -formula  CaCO3  1
    -m0       0                           
    -tol      1e-10
    -steps    600 in 100 steps
    -cvode    true

USER_GRAPH 1
-axis_titles "Time, sec" "Precipitate, mol"
10 GRAPH_X SIM_TIME/60
20 GRAPH_Y KIN("Calcite")
END

[dlparkhurst]

This script uses a threshold such that the surface-area-dependent rate doesn't kick in until there is a given amount of reactant present. It also allows for the possibility that at some point, calcite might dissolve again.

Code: [Select]

SOLUTION 1
    pH         7   charge
    density   1   calculate
    Cl         1
    Na         1
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Calcite   0   10
    CO2(g)     -1   
SAVE SOLUTION 1
END

USE SOLUTION 1

EQUILIBRIUM_PHASES
CO2(g) -3.4

INCREMENTAL_REACTIONS True

RATES
    Calcite
-start
 10 dpmet = ((1/TK)-(1/298.15))/8.3144
 20 srmin = SR("Calcite")
 30 moles = 0
 40 IF(M <= 0) AND (srmin < 1) THEN GOTO 180
 50 s = 0.7
 60 mm = 100.1
 80 knu = 0.000000148*EXP((-63500)*dpmet)
 90 k1 = 0.00658*EXP((-56100)*dpmet)*(ACT("HCO3-")^1.982)
100 kpre = (knu+k1)
110 theta = 0.5
120 eta = 2
130 threshold = 1e-9
140 IF(M <= 0 AND srmin <= 1) THEN GOTO 200
150 IF(M >= threshold) THEN GOTO 180
160 rate = 1e-9
170 GOTO 190
180 rate = s*M*mm*kpre*(ABS(1-srmin^theta)^eta)
190 moles = rate*TIME*(1-srmin)
200 SAVE moles
-end
   
KINETICS 1
Calcite
    -formula  CaCO3  1
    -m0       0                           
    -tol      1e-10
    -steps    600 in 100 steps
    -cvode    true

USER_GRAPH 1
    -headings               Time KIN(Calcite) SI(Calcite)
    -axis_titles            "Time, minutes" "Precipitate, mol" "SI(Calcite)"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X SIM_TIME/60
20 GRAPH_Y KIN("Calcite")
30 GRAPH_SY SI("Calcite")
  -end
    -active                 true
END

[MichaelZ20]

Thank you, David!