Calcite Dissolution Kinetics in PHREEQC Manual

[bschroth]

Hello:
I am running a model using the calcite dissolution RATE and KINETIC blocks from the Description of Data Input section of the PHREEQC manual, based on the cited reference of Plummer et al. (1978).  The water I am using is only mildly acidic (pH 6), but the dissolution rate I am getting is extremely low (order of 10^-10 mol/yr).  I adjusted the pH to lower values, and even at pH 3 the rate was still very low (order of 10^-7 mol/yr).  The example in the manual appears to be different from the Plummer et al. expression - is there an error in the manual's example, or is there an alternative expression that should be used?  I can send my input file for reference.  Thank you

Brian

[dlparkhurst]

The handling of surface area is different in the RATE definition for version 3 than for version 2. Look at the RATES definition in the phreeqc.dat file that you are using to see how the parameters are defined. I believe in version 3, surface area was changed to be area per mole; this allowed KINETICS blocks to be scaled more easily for PHAST.

[bschroth]

I see that the phreeqc.dat file has different values for the parameters.  I copied both the RATES and KINETICS blocks from that file into my model.  The rate definitely goes up; however, it is still much lower than expected.  For a pH 6 solution with SI_calcite at -2.3, a starting calcite mass of 4.5 moles loses 1.2 mmol over a 10-year period, with the solution kept at constant pH and CO2.  For the same solution titrated to pH 3 with HCl, there is a loss of 0.3 mol calcite over the same period.  Does this sound correct?

[dlparkhurst]

No, it does not sound correct. This is how I interpreted your conditions. You can argue about the surface area.

Code: [Select]

RATES
    Calcite
-start
  1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
  2   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 )
 80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME
200 SAVE moles
-end
END
PHASES
Fix_H+
    H+ = H+
    log_k     0
END
SOLUTION 1
    temp      25
    pH        6
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    C         1 Calcite    -2.3
    Ca        1
    -water    1 # kg
END
KINETICS 1
Calcite
    -formula  Calcite  1
    -m        4.5
    -m0       4.5
    -parms    100 1
    -tol      1e-08
-steps       86400 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
USE solution 1
EQUILIBRIUM_PHASES 1
    CO2(g)    -1.54 10
    Fix_H+    -6 HCl       10
USER_GRAPH 1
    -axis_titles            "Time, days" "Calcite dissolved, mmol" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/(3600*24)
20 GRAPH_y (4.5 - KIN("Calcite"))*1000
  -end
    -active                 true
END


[bschroth]

Your model produces more reasonable results.  The two main differences from my model were (1) integration time and (2) KINETICS block parameters.  For (1), you integrated over a 1-day period in 10 steps, while my time step was much longer - I think that was the main problem with my model.  What is the optimal time step to use for longer-duration simulations?  For (2), the two PARM values in your KINETICS block (100 m2 and 1.0) are quite different from those listed in the example block listed in phreeqc.dat (16.7 m2/mol and 0.6); I see what you mean about the area being up for interpretation.  100 m2/mol seems very high - for calcite is is 10,000 cm2/g, compared to 100 cm2/g assumed for the fine fraction of your 1978 study.  Is there a rule of thumb for choosing these parameters?

Using your settings, I changed the total time from 1 day to 30 days, with the same time step of 0.1 day.  My pH 6 water caused a dissolution of 0.15 moles in 30 days, and for comparison the pH 3 water dissolved nearly all 4.5 moles in 30 days.  Would this compare favorably to observed data?

Thanks for all your help on this!

[dlparkhurst]

The notes in RATES indicates parameter 1 is cm^2/mol.

In laboratories, calcite reacts to near equilibrium in hours to days, assuming sufficient surface area.

[bschroth]

Unfortunately I am still getting very slow loss of calcite.  I altered the problem to start with 0.028 mol of calcite, and ran for the same period and time step as your example, but after 1 day there was only about 0.1% dissolution.  The essential portions of my model code are listed below.  Where is my mistake?

PHASES
pH_fix
H+ = H+
logk  0.0

SOLUTION 1
temp   10   
pH   6.02   
pe   10   
redox   N(+5)/N(-3)   
units   mg/l   
density   1   
Na   5.93
K   4.27   
Ca   90.3   
Mg   15.1   
Cl   9.83   
Mn   6.39   
S(+6)   247   as SO4
Alkalinity   22.6   as Ca0.5(CO3)0.5
N(+5)   0.13   
N(-3)   0.045   
Si   11.1   as SiO2

PURE 1
CO2(g) -1.7 10.0
pH_fix -6.02 HCl 10.0


INCREMENTAL_REACTIONS true
RATES #use parameters from Plummer et al (1978) as shown in phreeqc.dat
   Calcite
   -start
1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
2   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 )
80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME # convert from mmol to mol
200 SAVE moles
   -end

KINETICS 1
   Calcite
     -formula CaCO3 0.9 MgCO3 0.1
     -m    0.028
     -m0   0.028
     -parms  100 1 
     -tol 1.e-8
-steps  86400 in 10 steps
-step_divide 1
-runge_kutta 3
-bad_step_max 500
END

[dlparkhurst]

I think it is all about surface area. In PWP 1978, pH-stat experiments used 20 g of calcite at ~50 cm^2/g. That is about 5000 cm^2/mol and 0.2 mol of calcite, or 1000 cm^2 per 0.8 L for the experiment. Your calculation has 2.8 cm^2 per liter (100x0.028). So your calculation is about 400 times slower.

[bschroth]

Using your values of 5000 cm2/mol for PARM(1), a PARM(2) of 1, and m0 of 0.2 mol calcite, I get approximately 6% loss of calcite (0.0126 mol) over the 24-hour simulation period, with pH fixed at 6.02 and log Pco2 fixed at -1.7.  I calculated the log dissolution rate (mmol/cm2/sec) to be about -6.8, and the set of curves in PWP shows perhaps -6.5 at this pH and Pco2 combination.  As it is in Region 3 of the graph, the rate would be expected to be slow; I don't think the "free drift" would apply with the pH fixed at 6.02.  Does this sound about right?

[dlparkhurst]

Seems pretty close to me. The actual surface area in the paper is about 44 cm^2/g, and are you comparing at the same temperature?

Free-drift allows the solution to approach equilibrium, so pH and solution composition change through the experiment.

[bschroth]

I made the adjustments for more precise SSA and temperature; the results were essentially the same.

The RELEASE file contains the following change for the 2015 version:
        --------
        svn 9462
        --------
        RATES definitions were revised in all databases. The area parameter
        was changed to specific area (area per mole of reactant). Also, for aqueous kinetic
        reactions, the solution volume was added explicitly. The changes
        were made to facilitate a new version of PHAST that uses a representative
        volume that does not have one liter of water.

How is the solution volume added explicitly?  If I have 16L of a solution to be added to a certain mass of calcite, how do I specify the solution volume?

[dlparkhurst]

There is a Basic function SOLN_VOL (with units liters) that is calculated when using the phreeqc.dat, Amm.dat, and pitzer.dat databases.

However, I don't think you need it for the calcite rate. That rate is dependent on the surface area and the saturation ratio, which incorporates the effects volume. You would use the solution volume if you had a rate expression with units of mol/L/s.

[pia said]

Hi, I read your post because I'm modeling kinetic dissolution of calcite in a batch experiment. I'm not clear with the values of parms. How i know the number that i have to choose?

[dlparkhurst]

The parameters for calcite differ between version 2 and version 3. You need to look at the documentation included in the RATES definition in the database file that you are using.

[pia said]

Thanks Mr. Parkhurst, I decided to use param(1) = 380000 because I have the area BET of my calcite = 0.38 m2/gr and param(2) = 0.6 because I looked in the database.

My code is

Code: [Select]

RATES
    Calcite
-start
  1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
  2   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 )
 80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME
200 SAVE moles
-end
END

SOLUTION 1
    temp      25
    pH        2
    pe        6.33 O2(g) -0.68
    redox     pe
    units     ppm
    density   1
    As(5)     4
    B         30
    Fe(2)     50
    S(6)      1000 charge
    -water    1 # kg
END

KINETICS 1
Calcite
    -formula  Calcite  1
    -m        0.0625
    -m0       0.0625
    -parms    380000 0.6
    -tol      1e-08
-steps        200 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
#INCREMENTAL_REACTIONS

USE solution 1

EQUILIBRIUM_PHASES 1
    CO2(g) -1.54 10
    Goethite 0 0   
    Fe(OH)3(a) 0 0
    #Calcite 0 10   

USER_GRAPH 1
    -headings               time CaCO3 pH
    -axis_titles            "Time, seg" "Calcite dissolved, mol" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_y (0.0625 - KIN("Calcite"))
#20 graph_y tot("Ca")
30 pH = -LA("H+")
40 graph_sy pH
  -end
    -active                 true

With the results i have a WARNING: negative moles in solution 1 for Ca. I don'k know how to fix this. Can someone help me please?

[dlparkhurst]

The Warnings do not mean that the calculation failed, only that some negative concentrations occurred during the integration. The integration method will automatically go back and use smaller time steps to avoid the negative concentrations.

Unless PHREEQC prints an "ERROR" message, the calculation should be correct.

[pia said]

thanks for your answer ! I have another question, m0 and m are moles or moles/L ?

[dlparkhurst]

There is a very nice description of the input for KINETICS in the manual for version 3.

Line 3: -m moles

moles --Current moles of reactant. As reactions occur, the moles will increase or decrease. Default is equal to initial moles if initial moles is defined, or 1.0 mol if initial moles is not defined. Optionally, m or -m .

Line 4: -m0 initial moles

initial moles --Initial moles of reactant. This identifier is useful if the rate of reaction is dependent on grain size. Formulations for this dependency often include the ratio of the amount of reactant remaining to the amount of reactant initially present. The quantity initial moles does not change as the kinetic reactions proceed. Frequently, the quantity initial moles is equal to moles at the beginning of a kinetic reaction. Default is equal to moles if moles is defined, or 1.0 mol if moles is not defined. Optionally, m0 or -m0

[pia said]

yes, I saw the manual but I don't know how phreeqc considers the volume of my solution

[dlparkhurst]

-water in SOLUTION

[pia said]

After doing the model in phreeqc I did an experimentation in laboratory, I dissolved 1.5 gr of calcite in acidic water (pH initial 2) and I measured the pH each 20 minutes. The results aren't consistent  whit my model because in the experiment the pH increment were slower than in the model, for example in the experiment the pH is almost 7 after 20 hours and in the model the pH is almost 7 after 3 hours. I’m frustrated because I thought the results would be similar but it didn’t happened. Now I’m trying to fix the model or explain why the results are not similar, if someone could help me I would be very grateful.

Code: [Select]

RATES
    Calcite
-start
  1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
  2   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 )
 80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
120 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O")) * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME
200 SAVE moles
-end

    Fe(OH)3(a)
-start
 10 REM parm(1)= A0 #m2/kg
 20 REM parm(2) = V #L
 30 rate = 10^-11 * (1 - SR("Fe(OH)3(a)")) * parm(1)/parm(2) * (m/m0)^0.67  #Grundl and Jim Delwiche
 40 moles = rate * time
200 SAVE moles
-end

    Gibbsite
-start
 10 REM parm(1)= A0 #m2/kg
 20 REM parm(2) = V #L
 30 rate = 10^-14 * (1 - SR("Gibbsite")) * parm(1)/parm(2) * (m/m0)^0.67  #verificar buscar referencia
 40 moles = rate * time
200 SAVE moles
-end

PHASES
Ca3(AsO4)2:4w       490
        Ca3(AsO4)2:4H2O = 3Ca+2 + 2AsO4-3 + 4H2O
        log_k           -18.905

Ca3(AsO4)2:3w      #(referencia de Ksp arsenatos de calcio)
        Ca3(AsO4)2:3H2O = 3Ca+2 + 2AsO4-3 + 3H2O
        log_k           -21.14

Ca5(AsO4)3OH
        Ca5(AsO4)3OH = 5Ca+2 + 3AsO4-3 + OH-
  log_k           -40.12

Ca4(OH)2(AsO4)2:4w       
        Ca4(OH)2(AsO4)2:4H2O = 4Ca+2 + 2AsO4-3 + 4H2O + 2OH-
        log_k           -27.49

SOLUTION 1  #acid water recipe
    units     ppm
    pH        2.0
    pe        6.33 O2(g) -0.68
    density   1
    Al        60.02
    As        2
    B         4
    Ca        265
    Cl        957
    Fe(3)     60.2
    K         70.3
    Mg        125
    Mn        10
    Na        351
    Pb        1
    S(6)      1883
    Zn        10
    -water    0.24 # kg se considera el volumen utilizado por daniela como disolucion de calita representativa

EQUILIBRIUM_PHASES
    CO2(g) -1.54 10
    O2(g) -0.68
   

KINETICS 1
Calcite
    -formula  Calcite  1
    -m0       0.015 #1.5 gr of calcite in solution
    -m        0.015
    -parms    380000 0.6 #area BET 0,38 m2/gr
    -tol      1e-08
Fe(OH)3(a)
    -formula  Fe(OH)3(a)  1 FeAsO4  0.052
    -parms     0.012 0.24
    -tol       1e-08
#Gibbsite
#    -formula  Gibbsite  1
    #-m        0.0005359
    #-m0       0.0005359
    #-parms    40 0.24
    #-tol      1e-08
-steps         36000 in 100 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
INCREMENTAL_REACTIONS True


USER_GRAPH 1
    -headings               Tiempo Ca pH
    -axis_titles            "Tiempo, horas" "Ca en solución, mg/l" ""
    -chart_title            "Modelo PWP sistema abierto, Datos Henry 2015"
    -axis_scale sy_axis     0 8 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_y TOT("Ca")*40.078*1000
30 pH = -LA("H+")
40 graph_sy pH
  -end
    -active                 true

USER_GRAPH 2
    -headings               _time_ Fe(3) Al pH
    -axis_titles            "Tiempo, horas" "Fe, mg/l" "pH" "Al, mg/l"
    -chart_title            "Concentración de Al y Fe en el timepo"
    -initial_solutions      true
    -axis_scale sy_axis     0 8
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y TOT("Fe(3)")*55.8*1000  TOT("Al")*27*1000
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

USER_GRAPH 3
    -headings               _time_ concentracion
    -axis_titles            "Tiempo, horas" "As mg/l" ""
    -chart_title            "As en el tiempo"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y TOT("As(5)")*75*1000
  -end
    -active                 true

USER_GRAPH 4
    -headings               _time_  Fe(3) #pH
    -axis_titles            "Tiempo, horas" "% remoción Fe disuelto" #"pH"
    -chart_title            "Remoción de Fe en función del tiempo"
    -initial_solutions      true
    -axis_scale sy_axis     0 8
    -axis_scale y_axis     0 120
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y (1 - TOT("Fe(3)")*55.8*1000/60.2)*100 #cambiar el hierro incial cuando sea necesario
#30 GRAPH_SY -LA("H+")
  -end
    -active                 true

USER_GRAPH 5
    -headings               _time_ As
    -axis_titles            "Tiempo, horas" "% remoción As" ""
    -chart_title            "Remoción de As en función del tiempo"
    -axis_scale y_axis      0 120 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y (1 - TOT("As(5)")*75*1000/2)*100
  -end
    -active                 true
END

[dlparkhurst]

It is not surprising that the kinetics calculations do not match your experiment exactly.

First, the rates of reaction are directly dependent on the surface area of calcite. Often surface area is used as a fit parameter to match the experimental data.

Another factor is mixing of the solution. Was the experiment stirred or not? resulting in a well-mixed solution or potentially quasi-diffusive transfer from the surface to the solution.

Other factors that can speed or retard the reaction are the rate of CO2 transfer to or from the solution, the temperature, purity of the calcite, any additional ions in solution, and more factors I haven't though of.

[pia said]

Thank you ver much !!

[pia said]

Hi again

When you said "Often surface area is used as a fit parameter to match the experimental data" how can I fit this parametrer; is manual ? or there is a preocedure ?

Thank you very much!

[dlparkhurst]

If you are just adjusting one parameter, you could do it manually. If you want to be a little more careful, you can use a fitting program. PHREEPLOT (http://phreeplot.org/), which is tightly linked to PHREEQC, could be used, or you can use a general parameter estimation program like PEST.