Error Reproducing "Fractionation of Carbon during Precipitation"

[GuyC]

Hi,

I am trying to reproduce the output from the example on p.592 from "Geochemistry, Groundwater and Pollution 2nd ed." under Numerical Modeling.

I am running through Mac OS X and am using the Phreeqc.dat database. I copied the example verbatim as shown below:

# Degas CaHCO3 solution, precipitate calcite, find d13C...

SOLUTION 1
 Ca 5 charge; C(4) 10


RATES
CO2(g)
-start
10 P_aq = sr("CO2(g)")
20 dCO2 = parm(2) * (10^parm(1) - P_aq) * time
30 save dCO2
40 put(dCO2, 5)
-end

Calcite

210 put(moles, 4) # store amount of calcite precipitate
-end

C13
# Distribute 13C among solute species 1 􏰘 CO2aq, 2 􏰘 HCO3-, 3 􏰘 CO3-2
# and in calcite and CO2 gas.
-start
10 aa_1_2 = exp-(9.866/TK - 24.12e-3)
20 aa_3_2 = exp(-0.867/TK + 2.52e-3)
30 aa_cc_2 = exp(-4.232/TK + 15.1e-3)
40 aa_g_2 = exp-(9.552/TK - 24.1e-3)

50 aH = act("H+")
60 K1 = aH * act("HCO3-") / act("CO2") / act("H2O")
70 K2 = aH * act("CO3-2") / act("HCO3-")
80 m13C = m/(aa_1_2*aH/K1+1+aa_3_2*K2/aH)
90 m12C = (tot("C(4)")) / (aH/K1 + 1 + K2/aH)

100 R2 = m13C / m12C
110 d13C = aa_cc_2 * -get(4) * R2
120 d13C_g = aa_g_2 * -get(5) * R2
130 save d13C + d13C_g
140 put(R2, 2); 180 put(d13C, 10); 190 put(d13C_g, 20)
-end

C13_cc; -start; 10 save -get(10); -end
C13_g; -start; 10 save -get(20); -end

KINETICS 1
 CO2(g);
        -m0 0;
        -parm -3.5 0.001
 Calcite;
        -m0 0;
        -parm 60 0.67

 #d13C = -10
 C13;
        -formula C 0;
        -m0 0.1112463e-3
 C13_cc;
        -formula C 0;
        -m0 0
 C13_g;
        -formula C 0;
        -m0 0
 -steps 1e-5 3 5 10 15 20 30 45 60 100 150 300 600 1e3 5e3 #time steps
INCREMENTAL_REACTIONS
END


Yet I get several iterations of the following output error:

WARNING: Negative moles in solution for C, nan. Recovering...

I've tried decreasing the target pCO2 parameter, removing the isotopes altogether, and changing the solution composition, but still cannot get rid of the error. The only thing I did not include was the USER_GRAPH section as I am running on Mac, but I don't think it contains any necessary functions to make a sensible output. What am I missing?

[dlparkhurst]

There is a difference between warnings and errors. Warnings do not cause the program to stop; in this case, the warning is saying that the kinetic integration is having a hard time, but it is trying to recover. The program can still produce an accurate results, even though it writes these warning messages.

Errors cause the program to stop. In general, the program should stop if it cannot produce an accurate answer for kinetics, or it cannot solve the nonlinear equations describing the chemical system.