Help with Kinetic reactions

[spence]

Hello all,

I've been working on a script to model my very first calcite dissolution by HCl acid. When I run the script with Phreeqc, my output data shows that concentration of H+ is zero for every time step, and my M calcite remained the same all through (meaning no reaction occurred). This seems a bit strange to me considering that I specified 0.01 M initial H+. To confirm what components I have in the system, I ran the same script with RM, and the GetComponents() function returns H, O, cb, C (which I suppose is HCO3-), Ca, and Cl, but no H+.

I assumed there is only pure water in the column prior to the influx of HCl. For my KINETICS, RATES, and TRANSPORT data blocks, I tried to follow the examples provided in the documentation. I suspect I may be doing something wrong, but I’ve tried and could not deduce what the issue might be – as I am new and still learning Phreeqc.

My input data is:

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.0-15749\database\pitzer.dat
TITLE Calcite dissolution attempt.
SOLUTION_MASTER_SPECIES
H H+ -1.0 H 1.008
H(0) H2 0.0 H
H(1) H+ -1.0 0.0
E e- 0.0 0.0 0.0
O H2O 0.0 O 16.0
O(0) O2 0.0 O
O(-2) H2O 0.0 0.0
C CO3-2 2 HCO3 12.0111
C(4) CO3-2 2 HCO3 12.0111

##############################
SOLUTION_SPECIES

H+ = H+
log_k 0.0
-gamma 9.0 0.0

e- = e-
log_k 0.0

H2O = H2O
log_k 0.0

H2O = OH- + H+
log_k -14.0
delta_h 13.362 kcal
-analytic -283.971 -0.05069842 13323.0 102.24447 -1119669.0
-gamma 3.5 0.0

2 H2O = O2 + 4 H+ + 4 e-
log_k -86.08
delta_h 134.79 kcal

2 H+ + 2 e- = H2
log_k -3.15
delta_h -1.759
           
##############################
SOLUTION 0  HCl
        units            mol/kgw
        temp             25.0
        pH               2.0     charge
        H(1)             0.01 # 0.01M HCl solution
        Cl               0.01
SOLUTION 1-4  Initial solution for column
        units            mol/kgw
        temp             25.0
        pH               7.0     charge
pure water       1.0
END
KINETICS 1-4
        Calcite
formula CaCO3 1
m 3.1
m0 3.1
parms 5.0 0.3
tol 1.e-8
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
TRANSPORT
        -cells           4
        -lengths         0.02
        -shifts          25
        -time_step       1440
        -flow_direction  forward
        -boundary_conditions   flux  flux
        -diffusion_coefficient 1.e-9
        -dispersivities  0.000
        -correct_disp    true
-punch_cells      4
        -punch_frequency  1
        -print_cells      4
        -print_frequency  1
USER_PUNCH
    -heading Calcite Act_H+
    -start
            10 punch M
    20 punch ACT("H(1)")
-end
SELECTED_OUTPUT
        -file            outputFile.sel
        -reset           false
        -step
        -totals          H(1) Cl Ca C
        -high_precision true
Thanks for any help, comments, or suggestions.

[Leo]

Hi,

I didn't check your whole model, but in your USER_PUNCH block, try using:

Code: [Select]

10 PUNCH KIN("Calcite")
20 PUNCH ACT("H+")

At least, this will give you the activity of H+ in the SELECTED_OUTPUT file. I tried this without a transport simulation and it works.

[dlparkhurst]

Here is a batch example of kinetic calcite dissolution as a function of time. I think the batch calculation is easier to understand without the added complexity of transport. When you move to transport, the calculation will reach a steady-state and not proceed past a certain point on the batch curve, depending on the length of the column and travel time.

Code: [Select]

#DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.0-15749\database\pitzer.dat
TITLE Calcite dissolution attempt.
         
##############################
SOLUTION 0  HCl
        units            mol/kgw
        temp             25.0
        pH               2.0     charge
        Cl               0.01
KINETICS 1
Calcite
formula CaCO3 1
m 3.1
m0 3.1
parms 5.0 0.3
tol 1.e-8
      -steps   86400 in 24
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
USE solution 0
USE kinetics 1
USER_GRAPH 1
    -headings               Time Delta_Calcite pH
    -axis_titles            "Time, hours" "Delta Calcite, moles" "pH"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  time
  -start
 5 GRAPH_X TOTAL_TIME/3600
10 GRAPH_Y KIN_DELTA("Calcite")
20 GRAPH_SY -LA("H+")
  -end


[spence]

Hi,

I didn't check your whole model, but in your USER_PUNCH block, try using:

Code: [Select]

10 PUNCH KIN("Calcite")
20 PUNCH ACT("H+")

At least, this will give you the activity of H+ in the SELECTED_OUTPUT file. I tried this without a transport simulation and it works.

Thank you @Leo for your suggestion. I appreciate.

[spence]

Thank you Dr Parkhurst. This will really help my understanding going forward. I really do appreciate.

I have a few questions thus far:
1. Is Delta_Calcite at every time step the difference between initial and the current moles of calcite? If true, I notice that when I manually subtract the current moles from the original, I get slightly different figures. I suppose this is due to round off errors?
2. How can I modify this script to simulate dissolution at equilibrium? I like to compare what happens if equilibrium (only) or kinetic (only) reactions occur.

Many thanks and regards.
Spence.

[dlparkhurst]

My bad. It is not roundoff error. The calculation was actually done twice, once because SOLUTION and KINETICS were defined before the END statement, and once because of the USE statements. Unlike other reactants, KINETICS is automatically updated and saved, so there was a decrease in the amount of calcite from the first calculation, and again for the second calculation. I had intended to add an END statement after the SOLUTION definition to perform the kinetic reaction only as a result of the USE statements.

For the pedantic, KIN_DELTA in this case gives the change in the kinetic reactant at each step relative to the initial state. There is an option INCREMENTAL_REACTIONS, in which case, KIN_DELTA would give the change in kinetic reactant between steps. INCREMENTAL_REACTIONS will speed up the calculation because it saves the results after each step and only integrates over the next time increment. Without INCREMENTAL_REACTIONS the integration is from the initial condition to the time of each step--beginning to time1, beginning to time2, etc.

The following dissolves calcite to equilibrium.

Code: [Select]

SOLUTION 0  HCl
        units            mol/kgw
        temp             25.0
        pH               2.0     charge
        Cl               0.01
END
EQUILIBRIUM_PHASES 1
Calcite 0 10
END

[spence]

Thank you Dr Parkhurst for your help thus far.

I have now setup the simulation using the TRANSPORT keyword (see attachment).

Considering the first paragraph in your last comment, I see that Delta Calcite = M – M0 for the batch reaction example. However, this equality does not hold for the transport case. Any possible explanation? I’m hoping there isn’t any double calculation happening in this case?

Also, my pH values are decreasing (see attached image) in contrast to the batch reaction example. My guess (ludicrous, I imagine) is that, in the case of the batch reaction, the outputted pH values are wrt H+ from the acid, which is being diluted as its concentration decreases, whereas in my transport case, pH is being computed wrt H+ from H20 (input solution), which is becoming acidic as the acid comes into the solution. Could this be right? If so, how do I extract the pH wrt H+ from the acid?

Many thanks and regards.
Spence.

[dlparkhurst]

No Dr please, I do not have a doctorate.

For TRANSPORT, KIN_DELTA("Calcite") returns the amount of calcite reacted over the transport step. The total amount of calcite reacted over the course of the simulation for a cell will be M0 - KIN("Calcite") assuming M0 = M initially (not necessarily the case). If you sum KIN_DELTA for each time step, it should be equal to M0 - KIN("Calcite").

The pH is slowly decreasing in the TRANSPORT case because the rate depends on M/M0, and M is slowly decreasing. Thus, the rate of reaction is slowly decreasing in each cell, and less reaction causes lower pH. If you had a very large amount for M0, the decrease would not be detectable at all.

The M/M0 factor is in the batch reaction also, but again the effect is small if M0 is large.

After one pore volume (4 time steps), your transport calculation is at steady state except for the small effect of M/M0.