Phreeqc seems to halve first kinetic step

[peterwadeuk]

Greetings good persons.

In the transport-kinetic script below I attempt to introduce a solution in steps of 30 days, to a 2-cell column.

One of my difficulties is that I seem to be receiving output from Phreeqc at periods of 15, 45, 75, 105 days. The start time is not zero or 30 days.

Grateful for any advice.

All the best,
Peter

========================================

TITLE Rebuild kinetic model of Mike's experiments
PHASES
Plagioclase
    Ca0.5Na0.5Al1.5Si2.5O8 + 6H+ + 2H2O = 1.5Al+3 + 0.5Ca+2 + 2.5H4SiO4 + 0.5Na+
    log_k     914.2 #     log_k     14.2 # Change to ensure reactivity
    delta_h   -12.467 kJ
Diopside
    CaMgSi2O6 + 4H+ + 2H2O = Ca+2 + 2H4SiO4 + Mg+2
    log_k     919.894 #    log_k     19.894 # Change to ensure reactivity
    delta_h   -32.348 kcal

SOLUTION 0 Pore waters in absence of basalt
    temp      25
    pH        8.156
    pe        12.471
    redox     pe
    units     mol/l
    density   1
Al   3.251E-07
Ba   5.471E-06
C   9.689E-04
Ca   7.100E-04
Cl   1.630E-04
F   5.841E-05
Fe   2.560E-08
K   2.010E-05
Li   3.360E-07
Mg   6.131E-05
Mn   4.011E-08
N   4.647E-04
Na   1.800E-04
S   5.081E-05
Si   5.389E-05
Sr   9.761E-07
Zn   2.310E-06
    -water    1 # kg

SOLUTION 1-2 Pore waters in absence of basalt
    temp      25
    pH        8.156
    pe        12.471
    redox     pe
    units     mol/l
    density   1
Al   3.251E-07
Ba   5.471E-06
C   9.689E-04
Ca   7.100E-04
Cl   1.630E-04
F   5.841E-05
Fe   2.560E-08
K   2.010E-05
Li   0
Mg   6.131E-05
Mn   4.011E-08
N   4.647E-04
Na   1.800E-04
S   5.081E-05
Si   5.389E-05
Sr   9.761E-07
Zn   2.310E-06
    -water    1 # kg

RATES
   Plagioclase
-start
10 moles = 0
20 IF(M <= 0 OR SI("Plagioclase")) >= 0 THEN GOTO 220
25 rock_area = PARM(1)
30 area_fract = PARM(2)
35 a0 = rock_area * area_fract
40 v = PARM(3)
50 dif_temp = 1/TK-1/298
60 k_acid = 10^-7.87
70 k_neut = 10^-10.91
80 k_base = 10^-15.57
90 eapp_acid = 42.1
100 eapp_neut = 45.2
110 eapp_base = 71
120 n_acid = 0.63
130 n_base = -0.57
150 hplus = MOL("H+")
160 r_acid = k_acid*EXP((-eapp_acid/8.314e-3)*dif_temp)*(hplus^n_acid)
170 r_neut = k_neut*EXP((-eapp_neut/8.314e-3)*dif_temp)*(hplus^n_neut)
180 r_base = k_base*EXP((-eapp_base/8.314e-3)*dif_temp)*(hplus^n_base)
190 r_all = r_acid+r_neut+r_base
200 rate = (a0/v)*(M/M0)^0.67*r_all*(1-SR("Plagioclase"))
205 print "r_acid, r_neut, r_base, SR-plag, rate, time = ", r_acid, r_neut, r_base, SR("Plagioclase"), rate, time
210 moles = rate*TIME
220 SAVE moles
-end
   Diopside
-start
10 moles = 0
20 IF(M <= 0 OR SI("Diopside")) >= 0 THEN GOTO 220
25 rock_area = PARM(1)
30 area_fract = PARM(2)
35 a0 = rock_area * area_fract
40 v = PARM(3)
50 dif_temp = 1/TK-1/298
60 k_acid = 10^-6.36
70 k_neut = 10^-11.11
80 k_base = 0
90 eapp_acid = 96.1
100 eapp_neut = 40.6
110 eapp_base = 0
120 n_acid = 0.71
130 n_base = 0
150 hplus = MOL("H+")
160 r_acid = k_acid*EXP((-eapp_acid/8.314e-3)*dif_temp)*(hplus^n_acid)
170 r_neut = k_neut*EXP((-eapp_neut/8.314e-3)*dif_temp)*(hplus^n_neut)
180 r_base = k_base*EXP((-eapp_base/8.314e-3)*dif_temp)*(hplus^n_base)
190 r_all = r_acid+r_neut+r_base
200 rate = (a0/v)*(M/M0)^0.67*r_all*(1-SR("Diopside"))
210 moles = rate*TIME
220 SAVE moles
-end

KINETICS 1
Plagioclase
    -formula  Ca0.5Na0.5Al1.5Si2.5O8Cs0.001  1
    -m        0.2789
    -m0       0.2789
    -parms    100.9 0.4661 1.131
    -tol      1e-08
Diopside
    -formula  CaMgSi2O6Rb0.001  1
    -m        0.0621
    -m0       0.0621
    -parms    100.9 0.0783 1.131
    -tol      1e-08
-steps       1
-step_divide 1
-runge_kutta 3
-bad_step_max 500

SELECTED_OUTPUT 1
    -file                 180419_1333_MikeExpts1_kinetic.tsv
    -totals               Na  Ca  Mg  C  Li Cs Rb
    -equilibrium_phases   Calcite  Kaolinite  CO2(g)
    -kinetic_reactants    Plagioclase Diopside
TRANSPORT
    -cells                 2
    -shifts                4
    -time_step             2592000 # seconds
    -lengths               2*0.25
    -diffusion_coefficient 0
    -thermal_diffusion     2   0
    -print_cells           1-2
    -punch_cells           1-2
EQUILIBRIUM_PHASES 1
    Calcite   0 10
    Kaolinite 3 0
    CO2(g)    -2.39 10
EQUILIBRIUM_PHASES 2
    Calcite   0 0
    Kaolinite 3 0
    CO2(g)    -3.39 10

[dlparkhurst]

I believe the calculation is accounting for the location of the node. The nodes are at half the length of the cell, so it is accounting for when a concentration front would reach the node.

[peterwadeuk]

Thank you David