Problems of Kinectics of calcite dissolution during reactive transport modelling

[Yongqiang]

Dear phreeqc experts,

I tried to use the following code to model reactive transport flow with the kinetics of calcite dissolution. However, my following code seems doesn't work. Can anyone help me out? Thank you very much.

Cheers,
Yongqiang

Code: [Select]

DATABASE c:\phreeqc\database\PITZER.DAT
TITLE reactive transport in carbonate reservoirs

SOLUTION 0
units mol/kgw
temp  25
pH 7
Na 3
Ca 0.05
Cl 3.06
# Na 3e-2
# Ca 3e-5
# Cl 3.06e-2
EQUILIBRIUM_PHASES 0
CO2(g) 2 100
#END

SOLUTION 1-20 formation brine
units mol/kgw
temp 25
pH 7
Na 3
Ca 0.03
Cl 3.06

EQUILIBRIUM_PHASES 1-20
Calcite
SAVE SOLUTION 1-20
END


RATES
Calcite
-start
1 rem M = current number of moles of calcite
2 rem M0 = number of moles of calcite initially present
3 rem PARM(1) = A/V, cm^2/L
4 rem 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 / TK )
40  k2 = 10^(2.84 - 2177.0 / TK)
50 if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
60 if TC > 25 then k3 = 10^(-1.1 - 1737.0 / TK )
70  t = 1
80 if M0 > 0 then t = M/M0
90 if t = 0 then t = 1
100 area = PARM(1) * (t)^PARM(2)
110 rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120 rem 1e-3 converts mmol to mol
130 rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
140 moles = rate * TIME
200 SAVE moles
-end

KINETICS 1-20
Calcite
-m 0
-m0 7.e-0
-parms 5.0 0.3
-tol 1.e-8

TRANSPORT
-cells 20
-lengths 2
-shifts 1
-time_step 0.020 minute
-flow_direction forward
-boundary_conditions flux flux
-diffusion_coefficient 0
-dispersivities 0
-correct_disp true
-punch_cells 1-20
-punch_frequency 1
# -print_cells 100
# -print_frequency 10

SELECTED_OUTPUT
-file carbonate reservoirs
-reset false
-step
-totals Na Cl K Ca
USER_PUNCH
-heading Ca
10 PUNCH TOT("Ca")
#20 PUNCH TOT("Calcite")


USER_GRAPH 1
-chart_title "CaCl Reactive Transport"
-headings pH Ca
-axis_titles "Distance" "Chemical species"
-axis_scale x_axis auto
-axis_scale y_axis auto
-initial_solutions false
-plot_concentration_vs x
-start
10 x = DIST
20 PLOT_XY x, -LA("H+")
30 PLOT_XY x, EQUI("Calcite")#SI("Calcite") #
#40 PLOT_XY x, TOT("Ca")
-end
END



[dlparkhurst]

Code: [Select]

TITLE reactive transport in carbonate reservoirs

RATES
Calcite
-start
1 rem M = current number of moles of calcite
2 rem M0 = number of moles of calcite initially present
3 rem PARM(1) = A/V, cm^2/L
4 rem 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 / TK )
40  k2 = 10^(2.84 - 2177.0 / TK)
50 if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
60 if TC > 25 then k3 = 10^(-1.1 - 1737.0 / TK )
70  t = 1
80 if M0 > 0 then t = M/M0
90 if t = 0 then t = 1
100 area = PARM(1) * (t)^PARM(2)
110 rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
120 rem 1e-3 converts mmol to mol
130 rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
140 moles = rate * TIME
200 SAVE moles
-end
END

SOLUTION 0
units mol/kgw
temp  25
pH 7
Na 3
Ca 0.05
Cl 3.06
EQUILIBRIUM_PHASES 0
CO2(g) 2 100
SAVE solution 0
END

SOLUTION 1 formation brine
units mol/kgw
temp 25
pH 7
Na 3
Ca 0.03
Cl 3.06
EQUILIBRIUM_PHASES 100
Calcite
SAVE SOLUTION 1-20
END

KINETICS 1-20
Calcite
-m 1
      -m0 7
-parms 5.0 1
END

TRANSPORT
-cells 20
-lengths 2
-shifts 20
-time_step 1 days
-flow_direction forward
-boundary_conditions flux flux
-dispersivities 0
-punch_frequency 20
      -print_frequency 20

USER_GRAPH 1
-headings Dist pH SI("Calcite")
-axis_titles "Distance" "pH" "SI"
-start
10 GRAPH_X DIST
20 GRAPH_Y -LA("H+")
30 GRAPH_SY SI("Calcite")
-end
END


[Yongqiang]

Professor Parkhurst, thank you, you are really the king of hydrogeology!!

Thank you,
Yongqiang