equilibrium or kinetics

[bessie000]

I want to set up a model, which is composed of several minerals. They would first equilibrate with rain water. Then there is a boundary water, flowing through the assemblage, considering kinetic reactions among most of them. However, there were few minerals completely dissolved during the equilibrium section (magnetite, amorphous, hannebachite, ettringite), even though they were redefined in the kinetics block. Can some take a look and see where I got wrong?
Here is what I have:
TITLE Simulation 1.--Initial minerals equilibrate with sythetic groundwater

PHASES
Ettringite
Ca6Al2(SO4)3(OH)12:26H2O + 12 H+ = 6 Ca+2 + 2 Al+3 + 3SO4-2 + 38 H2O; log_k 56.4

Amorphous_phase
Si0.761Al0.15Ca0.0012K0.004Mg0.0011Na0.0076O1.5301(OH)0.45 + 0.4662 H+ + 1.0639 H2O = 0.0012 Ca+2 + 0.004 K+ + 0.0011 Mg+2 + 0.0076Na+ + 0.761 H4SiO4 + 0.15 Al+3; log_k 0.06271

Magnetite
Fe3O4 + 8 H+ = 2 Fe+3 + Fe+2 + 4 H2O; log_k 3.737

Mullite
Al6Si2O13 + 18 H+ = 6 Al+3 + 5 H2O + 2 H4SiO4; log_k 55

Hannebachite
CaSO3:0.5H2O + 0.5 O2 = Ca+2 + SO4-2 + 0.5 H2O; log_k 38.5352

Portlandite
Ca(OH)2 + 2 H+ = Ca+2 + 2 H2O; log_k 22.5552


SOLUTION 1-100 concentrated rainwater from Example 14
         pH 4.6
         pe 4.0   O2(g) -0.7
         temp 25
         units mmol/l
         Ca  0.191625
         Mg  0.035797
         Na  0.122668
         Cl  0.133704
         C   0.1096
         S   0.235153  charge
   
EQUILIBRIUM_PHASES 1-100
 
    Calcite   0 0.015
    Gypsum    0 0.020  #treat as in equilibrium in simulation for now due to its fast rate k25
    Hematite  0 0.039
    Quartz    0 0.074
    Magnetite 0 0.009
    Mullite   0 0.029
    Amorphous_phase 0 0.178
    Hannebachite   0  0.384
    Ettringite     0  0.0007
    O2(g)     -0.67 10   #partial pressure of O2 in dry air is 0.209 atm
    CO2(g)    -3.5  10   #partial pressure of CO2 in dry air is 3.16E-4

 
SAVE Solution 1-100
END

TITLE Simualation 2.---Kinetic and Rate Law Blocks

SOLUTION 0 #define recharge water here, use synthetic groundwater composition

    temp      25
    pH        5.8
    pe        4
    redox     pe
    units     mol/l
    density   1     
    Ca        0.008025
    O(0)      1.63e-4
    S(6)      0.0093125
    Mg        0.00307
    K         1.278e-5
    Na        4.34e-5
    Fe(2)     1.428e-7
    Zn        3.33e-7
    Al        0.00111
    Mn(2)     1.1376e-4

KINETICS 1-100 # used to specify kinetic reactions and parameters for batch-reaction and reactive-transport calculations
Calcite
-m   0.015  #moles, current moles of reactant
-m0  0.015  #initial moles of reactant
-tol 1e-6

Hematite
-m 0.039
-m0 0.039
-tol 1e-6

Quartz
-m 0.074
-m0 0.074
-tol 1e-6

Magnetite
-m 0.009
-m0 0.009
-tol 1e-6

Amorphous_phase
-m 0.178
-m0 0.178
-tol 1e-6

Hannebachite
-m 0.384
-m0 0.384
-tol 1e-6

Ettringite
-m 0.0007
-m0 0.0007
-tol 1e-6

Mullite
-m 0.029
-m0 0.029
-tol 1e-4

RATES #used to define mathematical rate expressions for kinetic reactions, assume 25 C, therefore the temperature term goes to 0 in the rate equation
    Calcite
    -start
    10 rem use rate law from Steefe and Lasaga (1994)
    20 rem rate parameters from Palandri and Kharaka (2004), suface area A is 0.0458 m^2/g
    30 si_cal = SI("Calcite")              #define saturation index (log(Q/K))for rate equation
    40 Q_K = 10^si_cal                      #calculate Q/K
    50 k_nu = 1.54882e-6                    #mol/m2.s, kinetic rate constant at neutral condition
    60 k_h = 0.501187 * (ACT("H+"))^1         #kinetic rate constant at acid condition
    70 k_oh = 0.000331131 * (ACT("OH-"))^1    #kinetic rate constant at base condition
    80 k = k_nu + k_h + k_oh                #k is kinetic rate constant
    90 r = k * 0.0458 * (1-Q_K)             #r is kinetic rate in mol/g.s, equals to k*A*(1-Q/K)
    100 rem converts r into mol/s by multiplying mass of calcite in the system
    110 r = r * 1.5 #mass of calcite in current set up is 1.5 g
    120 moles = r * TIME
    130 SAVE moles
    -end
   
Hematite
   -start
   10 rem use rate law from Steefe and Lasaga (1994)
   20 rem rate parameters from Palandri and Kharaka (2004), surface area A is 0.00129 m^2/g
   30 si_hem = SI ("Hematite")    #define saturation index log(Q/K) for rate equation
   40 Q_K = 10^si_hem             #calculate Q/K
   50 k_nu = 2.512e-15                    #mol/m2.s, kinetic rate constant at neutral condition
   60 k_h = (4.074e-10) * (ACT("H+"))^1         #kinetic rate constant at acid condition
   70 k_oh = 0                               #kinetic rate constant at base condition
   80 k = k_nu + k_h + k_oh                #k is kinetic rate constant
   90 r = k * 0.00129 * (1-Q_K)             #r is kinetic rate in mol/g.s
   100 rem converts r into mol/s by multiplying mass of hematite in the system
   110 r = r * 6.25 #mass of calcite in current set up is 6.25 g
   120 moles = r * TIME
   130 SAVE moles
   -end

 Quartz
   -start
   10 rem use rate law from Steefe and Lasaga (1994)
   20 rem rate parameters from Palandri and Kharaka (2004), suface area A is 0.00091 m^2/g
   30 si_qtz = SI("Quartz")             #define saturation index log(Q/K) for rate equation
   40 Q_K = 10^si_qtz                   #Q_K is Q/K
   50 k_nu = 1.023e-14                    #mol/m2.s, kinetic rate constant at neutral condition,no acid or base mechanism for quartz
   60 k = k_nu                            #k is kinetic rate constant
   70 r = k * 0.0458 * (1-Q_K)             #r is kinetic rate in mol/g.s
   80 rem converts r into mol/s by multiplying mass of quartz in the system
   90 r = r * 4.4 #mass of calcite in current set up is 4.4 g
   100 moles = r * TIME
   110 SAVE moles
    -end
   
Magnetite
   -start
   10 rem use rate law from Steefe and Lasaga (1994)
   20 rem rate parameters from Palandri and Kharaka (2004), surface area A is 0.019 m^2/g, Mayant et al. 2008
   30 si_mag = SI("Magnetite")
   40 Q_K = 10^si_mag
   50 k_nu = 1.660e-11                 #mol/m2.s
   60 k_h = (4.467e-9) * (ACT("H+"))^0.279         #kinetic rate constant at acid condition
   70 k_oh = 0                               #kinetic rate constant at base condition
   80 k = k_nu + k_h + k_oh                #k is kinetic rate constant
   90 r = k * 0.019 * (1-Q_K)             #r is kinetic rate in mol/g.s
   100 rem converts r into mol/s by multiplying mass of magnetite in the system
   110 r = r * 2.075 #mass of calcite in current set up is 2.075 g
   120 moles = r * TIME
   130 SAVE moles
   -end

Amorphous_phase 
   -start
   10 rem use rate law from Steefe and Lasaga (1994)
   20 rem rate parameters from Palandri and Kharaka (2004), surface area A is 0.62 m^2/g, Oelkers and Gislason, 2001
   30 si_amor = SI ("Amorphous_phase")
   40 Q_K =  10^si_amor
   50 k_nu = 2.23872e-12                  #mol/m2.s, kinetic rate constant at neutral condition
   60 k = k_nu                            #k is kinetic rate constant, no acid or base mechanism
   70 r = k * 0.62 * (1-Q_K)             #r is kinetic rate in mol/g.s
   80 rem converts r into mol/s by multiplying mass of amorphous in the system
   90 r = r * 89.62 #mass of calcite in current set up is 89.62 g
   100 moles = r * TIME
   110 SAVE moles
    -end

Hannebachite
   -start
   10 rem use rate law from Steefe and Lasaga (1994)
   20 rem rate parameters from Palandri and Kharaka (2004), surface area A is 0.071 m^2/g, Taerakul et al., 2003.
   30 si_han = SI ("Hannebachite")
   40 Q_K =  10^si_han
   50 k_nu = 2.733e-4   #use CaSO3 as analog,Hao and Dick 2000
   60 k = k_nu          #k is kinetic rate constant, no acid or base mechanism
   70 r = k * 0.071 * (1-Q_K)
   80 rem converts r into mol/s by multiplying mass of hannebachite in the system
   90 r = r * 49.5 #mass of calcite in current set up is 89.62 g
   100 moles = r * TIME
   110 SAVE moles
    -end

Ettringite
   -start
   10 rem use rate law from Steefe and Lasaga (1994)
   20 rem rate parameters from Palandri and Kharaka (2004), surface area A is 0.098 m^2/g, Baur et al. 2014
   30 si_ett = SI("Ettringite")
   40 Q_K = 10^si_ett
   50 k_nu = 1.14e-12
   60 k = k_nu
   70 r = k * 0.098 * (1-Q_K)
   80 rem converts r into mol/s by multiplying mass of ettringite in the system
   90 r = r * 0.883 #mass of calcite in current set up is 0.883 g
   100 moles = r * TIME
   110 SAVE moles
   -end

Mullite
   -start
   10 rem use rate law from Steefe and Lasaga (1994)
   20 rem rate parameters from Palandri and Kharaka (2004), surface area A is 0.124 m^2/g
   30 si_mul = SI ("Mullite")
   40 Q_K = 10^si_mul                   #Q_K is Q/K
   50 k_nu = 1.023e-14                    #mol/m2.s, kinetic rate constant at neutral condition
   60 k = k_nu                            #k is kinetic rate constant
   70 r = k * 0.124 * (1-Q_K)             #r is kinetic rate in mol/g.s
   80 rem converts r into mol/s by multiplying mass of mullite in the system
   90 r = r * 12.36 #mass of mullite in current set up is 12.36 g
   100 moles = r * TIME
   110 SAVE moles
    -end         



PRINT
-selected_output true
-status false


TRANSPORT
-cells 100
-length 1.0
-shifts 200
-punch 10

USER_GRAPH
-plot_concentration_vs time
-axis_title "pore volume or shift number" "concentration in mol/l" "pH"
-start
10 graph_x step_no #step_no is shifts
20 graph_y tot("Ca") #red
30 graph_SY -LA("H+")  #green

SELECTED_OUTPUT 
        -file           Finalsolution.sel
        -reset false
        -step
        -equilibrium_phases  Calcite Gypsum Hematite Amorphous_phase Quartz Mullite Hannebachite Magnetite Ettringite
       
END

[dlparkhurst]

You have defined EQUILIBRIUM_PHASES 1-100 and KINETICS 1-100. Both will be active during the TRANSPORT calculation. Regardless of the kinetic reactions, the phases in equilibrium phases will ensure the resulting solution is in equilibrium with the minerals, provided some moles of the minerals remain to react.

If you want to switch to purely kinetic reactions, you should DELETE equilibrium_phases 1-100. You don't really need to define 100 equilibrium phases if you just want to establish an initial solution composition. You could define an equilibrium phases with a number that is not in the column (0 or > 100) and react it with a solution, then SAVE that reacted solution as 1-100.

[bessie000]

Thank you so much,  Mr Parkhurst! That explains why. I was thinking the reactions occur in the order I write them, so first equilibrium with rain water, then kinetic reactions with boundary water, I did not realize they would both be active in the transport section.

[bessie000]

You have defined EQUILIBRIUM_PHASES 1-100 and KINETICS 1-100. Both will be active during the TRANSPORT calculation. Regardless of the kinetic reactions, the phases in equilibrium phases will ensure the resulting solution is in equilibrium with the minerals, provided some moles of the minerals remain to react.

If you want to switch to purely kinetic reactions, you should DELETE equilibrium_phases 1-100. You don't really need to define 100 equilibrium phases if you just want to establish an initial solution composition. You could define an equilibrium phases with a number that is not in the column (0 or > 100) and react it with a solution, then SAVE that reacted solution as 1-100.

Mr Parkhurst, I have tried what you suggested above, by making a new equilibrium phase 0 as the mineral assemblage, reacting with a rainwater, then save the solution as 1-100 as initial solution. Then I did pure kinetic reactions 1-100 in the TRANSPORT block. However, when I check my output files, all the minerals become completely dissolved, and this is very strange.

[dlparkhurst]

I suggest you use one SOLUTION and one KINETICS reactions to define a batch reactor, and run the reactor for a series of times, say 1 10 100 1000 ... seconds with INCREMENTAL_REACTIONS.

It is much easier to try to determine the problem using a single cell.

You have CaSO3, which may be the problem. SO3 is unstable relative to S-2 and SO4-2 in most cases, and I suspect this will cause all of the Hannebachite will disappear.

[bessie000]

Thanks, I will try what you suggested above to define a batch-reactor. But just being very curious, if I define a time_step above, even though very small (1440s), it is very slow to calculate, it took a whole day to calculate the chemical changes simulated in a few days. I think this is not normal. I have tried to eliminate all kinetics, it runs much faster. Does this mean there is a problem in my code or I used too many kinetic reactions?

[dlparkhurst]

Try using -cvode in KINETICS. It invokes an implicit solver for the kinetics that can be faster for "stiff" kinetic reactions. Also use INCREMENTAL_REACTIONS to avoid reintegrating the preceding time periods.

[Martin K]

Hi,
I have a note not exactly regarding your question. You are using kinetic data from Palandri and Kharaka (2004) where the carbonate mechanism instead of base mechanism is introduced in the rate equation for carbonate minerals. So I think ACT(“OH“) should be replaced with ACT(“CO2“) in the expression for k_oh in your code. Or am I wrong?
This issue is evoking my next question for all – is the activity of CO2 the proper thermodynamic quantity to be used in the rate equation for carbonate minerals used by Palandri and Kharaka (A Compilations of Rate Parameters …, USGS, 2004, eq. 7)? Or should it be ACT(“HCO3-”), PR_P(“CO2(g)”) or something like that?
Thank you,
Martin

[dlparkhurst]

I never noticed the distinction before, but you are right. For carbonate minerals, the base mechanism is replaced with the "Carbonate Mechanism" in table 33.

It says "reaction order n with respect to P(CO2)", although it also says the "rate parameters were regressed using eqn 11", which has a pH dependence. Seeing as they specifically say "Carbonate Mechanism", I assume the act(H+) factor in the basic mechanism is replaced with P(CO2). Is that your take?

For PHREEQC, ignoring the difference between partial pressure and fugacity, the PCO2 is either SR("CO2(g)") or 10^SI("CO2(g").

[Martin K]

Assuming that P(CO2) is probably easier to control than activity of CO2(aq) or HCO3- when conducting rate experiments on "carbonate mechanism", I think P(CO2) is the right variable.
Thank you for your insight.