kinetic of calcite dissolution

[P.Merdy]

I try to simulate the dissolution of calcite over time. I have been struggling with the following script for hours without undersdanding what's wrong. I have an error message "Moles of reaction not SAVEed for Calcite_kinetic" even thougt I wrote "SAVE moles" in the script. Could you help me to understand my errors, please? Thank you very much in advance.

SOLUTION 1
temp      25   
pH        7   
pe        4   
units     mol/kgw   
Ca        0   
C(4)      0

EQUILIBRIUM_PHASES 1   
Calcite 0  1

KINETICS 1
Calcite_kinetic
-formula CaCO3 1   
-m0 1.0   
-parms 1e-5 1.0   
-steps 86400 in 10 steps # in seconds
-cvode true

RATES
Calcite_kinetic
-start 
10 k = PARM(1)
20 n = PARM(2) 
30 si_calcite = SI("Calcite") 
40 omega = 10^si_calcite 
50 if (omega <= 1) then
60 rate = 0 
70 else
80 rate=k*(1-1/omega)^n 
90 end if
100 moles = rate * TIME 
110 SAVE moles
-end

SELECTED_OUTPUT   
-file Calcite_kinetic.sel   
-reset false   
-time true   
-totals Ca C(4)
-user_punch true

USER_PUNCH
-headings Step Time SI_Calcite Omega_Calcite Moles_Calcite
10 PUNCH STEP_NO, TIME, SI("Calcite"), 10^SI("Calcite"), KIN("Calcite_kinetic")

USER_GRAPH 1
20 GRAPH_X time
30 GRAPH_Y TOT("Ca")
-end
END

[dlparkhurst]

The Basic interpreter is pretty basic. The syntax for IF statements is

IF ... THEN ... ELSE ...

You cannot continue on the next numbered line, and there is no END IF statement.

When the interpreter gets to the END IF, it interprets it as END, so it never gets to the SAVE line.

Code: [Select]

SOLUTION 1
temp      25   
pH        7   
pe        4   
units     mol/kgw   
Ca        0   
C(4)      0

EQUILIBRIUM_PHASES 1   
Calcite 0  1

KINETICS 1
Calcite_kinetic
-formula CaCO3 1   
-m0 1.0   
-parms 1e-5 1.0   
-steps 86400 in 10 steps # in seconds
-cvode true

RATES
Calcite_kinetic
-start
10 k = PARM(1)
20 n = PARM(2)
30 si_calcite = SI("Calcite")
40 omega = 10^si_calcite
50 if (omega <= 1) then rate = 0 else rate=k*(1-1/omega)^n
100 moles = rate * TIME
110 SAVE moles

-end

SELECTED_OUTPUT   
-file Calcite_kinetic.sel   
-reset false   
-time true   
-totals Ca C(4)
-user_punch true

USER_PUNCH
-headings Step Time SI_Calcite Omega_Calcite Moles_Calcite
10 PUNCH STEP_NO, TIME, SI("Calcite"), 10^SI("Calcite"), KIN("Calcite_kinetic")

USER_GRAPH 1
20 GRAPH_X time
30 GRAPH_Y TOT("Ca")
-end
END

[P.Merdy]

Thank you so much for the explanations and your prompt response!

[P.Merdy]

I applied some corrections to my script, by removing calcite from 'equilibrium_phase' and replacing by phases instead. I chose longer times to apply the dissolution calculus. I would like the graph dissolved calcite vs time. The graph showed 0 in x-axis all along the time whatever the time scale I choose. If I choose "moles" or "tot("Ca"), it does change anything (I got zero in x-axis) as if calcite does not dissolve. I don't see my mistakes. Could you help me?

SOLUTION 1
temp      25   
pH        7   
pe        4   
units     mol/kgw   
Ca        0   
C(4)      0

PHASES 1   
Calcite 0  1

KINETICS 1
Calcite_kinetic
-formula CaCO3 1   
-m0 1.0   
-parms 1e-5 1
-steps 864000000 in 10 steps # in seconds
-cvode true

RATES
Calcite_kinetic
-start
10 k = PARM(1)
20 n = PARM(2)
30 si_calcite = SI("Calcite")
40 omega = 10^si_calcite
50 if (omega <= 1) then rate = 0 else rate=k*(1-1/omega)^n
100 moles = rate * TIME
110 SAVE moles

-end

SELECTED_OUTPUT   
-file Calcite_kinetic.sel   
-reset false   
-time true   
-totals Ca C(4)
-user_punch true

USER_PUNCH
-headings Step Time SI_Calcite Omega_Calcite Moles_Calcite
10 PUNCH STEP_NO, TIME, SI("Calcite"), 10^SI("Calcite"), KIN("Calcite_kinetic")

USER_GRAPH 1
20 GRAPH_X total_time
30 GRAPH_Y moles
-end
END
 

[dlparkhurst]

Think about it:

Code: [Select]

50 if (omega <= 1) then rate = 0 else rate=k*(1-1/omega)^n
Also, moles is not defined in USER_GRAPH.

[P.Merdy]

I thought about it. I should consider the dissolution rate of calcite as following
50 if (omega => 1) then rate = 0 else rate=k*(1-omega)^n

In other words, I have dissolution if omega is < 1 and I changed the equation. Hopefully, it's correct now. 

[dlparkhurst]

Your solution was undersaturated, so there was no dissolution. The logic was backwards.

[P.Merdy]

Thank you very much for your help, I got it.