| Next|| Previous || Top |

Example 6.-- Reaction-Path Calculations

In this example, the precipitation of phases as a result of incongruent dissolution of K-feldspar (microcline) is investigated. Only a limited set of phases--K-feldspar, gibbsite, kaolinite, and K-mica (muscovite)--is considered in this example. The reaction path for this set of phases was originally addressed by Helgeson and others (1969). In this example, the thermodynamic data for the phases (table 21, PHASES keyword) are derived from Robie and others (1978) and are the same as test problem 5 in the PHREEQE manual (Parkhurst and others, 1980).

PHREEQC can be used to solve this problem in three ways: the individual intersections of the reaction path and the phase boundaries on a phase diagram can be calculated (example 6A), the reaction path can be calculated incrementally (6B), or the reaction path can be calculated as a kinetic process (6C). In the first approach, no knowledge of the amounts of reaction is needed, but a number of simulations are necessary to find the appropriate phase-boundary intersections. In the second approach, only one simulation is sufficient, but the appropriate amounts of reaction must be known beforehand. In the third approach, a kinetic rate expression is used to calculate the reaction path, using a step-size adjusting algorithm which takes care of phase boundary transitions by automatically decreasing the time interval when necessary. Only the total time to arrive at the point of K-feldspar equilibrium is required. All three approaches are demonstrated in this example. PHREEQC implicitly contains all the logic of a complete reaction-path program (for example Helgeson and others, 1970, Wolery, 1979, Wolery and others, 1990). Moreover, the capability to calculate directly the phase boundary intersections provides an efficient way to outline reaction paths on phase diagrams, and the option to add the reaction incrementally and automatically find the stable phase assemblage allows points on the reaction path between phase boundaries to be calculated easily and rapidly. The kinetic approach and the Basic interpreter that is embedded in PHREEQC can be used to save and print the arrival time and the aqueous composition at each phase transition.

Conceptually, the example considers the reactions that would occur if K-feldspar were placed in a beaker and allowed to react slowly. As K-feldspar dissolves, other phases may begin to precipitate. In this example, it is assumed that only gibbsite, kaolinite, or K-mica can form, and that these phases will precipitate reversibly if they reach saturation. Phases precipitated at the beginning of the reaction may redissolve as the reaction proceeds.

Table 21. --Input data set for example 6

TITLE Example 6A.--React to phase boundaries.
            pH    7.0 charge
          temp   25.0
                Al(OH)3 + 3 H+ = Al+3 + 3 H2O
                log_k          8.049
                delta_h        -22.792 kcal
                Al2Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 2 Al+3
                log_k          5.708
                delta_h        -35.306 kcal
                KAl3Si3O10(OH)2 + 10 H+ = 3 Al+3 + 3 H4SiO4 + K+
                log_k          12.970
                delta_h        -59.377 kcal
                KAlSi3O8 + 4 H2O + 4 H+ = Al+3 + 3 H4SiO4 + K+
                log_k          0.875
                delta_h        -12.467 kcal
        -file           ex6A-B.sel
        -activities     K+ H+ H4SiO4
        -si             Gibbsite Kaolinite K-mica K-feldspar
        -equilibrium    Gibbsite Kaolinite K-mica K-feldspar
TITLE Example 6A1.--Find amount of K-feldspar dissolved to 
                    reach gibbsite saturation.
USE solution 1
        Gibbsite        0.0     KAlSi3O8        10.0
        Kaolinite       0.0     0.0
        K-mica          0.0     0.0
        K-feldspar      0.0     0.0
TITLE Example 6A2.--Find amount of K-feldspar dissolved to 
                    reach kaolinite saturation.
USE solution 1
        Gibbsite        0.0     0.0
        Kaolinite       0.0     KAlSi3O8        10.0
        K-mica          0.0     0.0
        K-feldspar      0.0     0.0
TITLE Example 6A3.--Find amount of K-feldspar dissolved to 
                    reach K-mica saturation.
USE solution 1
        Gibbsite        0.0     0.0
        Kaolinite       0.0     0.0
        K-mica          0.0     KAlSi3O8        10.0
        K-feldspar      0.0     0.0
TITLE Example 6A4.--Find amount of K-feldspar dissolved to 
                    reach K-feldspar saturation.
USE solution 1
        Gibbsite        0.0     0.0
        Kaolinite       0.0     0.0
        K-mica          0.0     0.0
        K-feldspar      0.0     KAlSi3O8        10.0
TITLE Example 6A5.--Find point with kaolinite present, 
                    but no gibbsite.
USE solution 1
        Gibbsite        0.0     KAlSi3O8        10.0
        Kaolinite       0.0     1.0
TITLE Example 6A6.--Find point with K-mica present, 
                    but no kaolinite
USE solution 1
        Kaolinite       0.0     KAlSi3O8        10.0
        K-mica          0.0     1.0
TITLE Example 6B.--Path between phase boundaries.
USE solution 1
        Kaolinite       0.0     0.0
        Gibbsite        0.0     0.0
        K-mica          0.0     0.0
        K-feldspar      0.0     0.0
        K-feldspar      1.0
        0.04 0.08 0.16 0.32 0.64 1.0 2.0 4.0 
        8.0 16.0 32.0 64.0 100 200 umol
TITLE Example 6C.--kinetic calculation
        -units mol/kgw
        Al         1.e-13
        K          1.e-13
        Si         3.e-13
        Gibbsite   0.0  0.0
        Kaolinite  0.0  0.0
        K-mica     0.0  0.0
#        k0 * A/V = 1e-16 mol/cm2/s * (10% fsp, 0.1mm cubes) 136/cm = 136.e-13 mol/dm3/s
        -parms 1.36e-11
        -m0 2.16
        -m  1.94
        -step_divide 1e-6
        -steps    1e2 1e3 1e4 1e5 1e6 1e7 1e8
  10  REM store the initial amount of K-feldspar
  20  IF EXISTS(1) = 0 THEN PUT(M, 1)
  30  REM calculate moles of reaction
  40  SR_kfld = SR("K-feldspar")
  50  moles = PARM(1) * (M/M0)^0.67 * (1 - SR_kfld) * TIME
  60  REM The following is for printout of phase transitions
  80  REM      Start Gibbsite
  90  if ABS(SI("Gibbsite")) > 1e-3 THEN GOTO 150
  100   i = 2
  110   GOSUB 1500
  150 REM      Start Gibbsite -> Kaolinite
  160 if ABS(SI("Kaolinite")) > 1e-3 THEN GOTO 200
  170   i = 3
  180   GOSUB 1500
  200 REM      End Gibbsite -> Kaolinite
  210 if ABS(SI("Kaolinite")) > 1e-3 OR EQUI("Gibbsite") > 0 THEN GOTO 250
  220   i = 4
  230   GOSUB 1500
  250 REM      Start Kaolinite -> K-mica
  260 if ABS(SI("K-mica")) > 1e-3 THEN GOTO 300
  270   i = 5
  280   GOSUB 1500
  300 REM      End Kaolinite -> K-mica
  310 if ABS(SI("K-mica")) > 1e-3 OR EQUI("Kaolinite") > 0 THEN GOTO 350
  320   i = 6
  330   GOSUB 1500
  350 REM      Start K-mica -> K-feldspar
  360 if ABS(SI("K-feldspar")) > 1e-3 THEN GOTO 1000
  370   i = 7
  380   GOSUB 1500
  1000 SAVE moles
  1010 END
  1500 REM subroutine to store data
  1510 if GET(i) >= M THEN RETURN
  1520     PUT(M, i)
  1530     PUT(TOTAL_TIME, i, 1)
  1540     PUT(LA("K+")-LA("H+"), i, 2)
  1550     PUT(LA("H4SiO4"), i, 3)
  1560 RETURN  
  10 DATA "A: Gibbsite               ", "B: Gibbsite  -> Kaolinite ", \
          "C: Gibbsite  -> Kaolinite ", "D: Kaolinite -> K-mica    ", \
          "E: Kaolinite -> K-mica    ", "F: K-mica    -> K-feldspar"
  20 PRINT "        Transition                 Time   K-feldspar      LA(K/H)   LA(H4SiO4)"
  30 PRINT "                                             reacted"
  40 PRINT "                                             (moles)"
  50 FOR i = 2 TO 7
  60   READ s$
  70   IF EXISTS(i) THEN PRINT s$, GET(i,1), GET(1) - GET(i), GET(i,2), GET(i,3)
  80 NEXT i
        -file ex6C.sel
        -reset false
  -head pH+log[K]  log[H4SiO4]
  10 PUNCH LA("K+")-LA("H+") LA("H4SiO4")

The input data set (table 21) first defines pure water with SOLUTION input and the thermodynamic data for the phases with PHASES input. Some of the minerals are defined in the database file ( phreeqc.dat ), but inclusion in the input data set replaces any previous definitions for the duration of the run (the database file is not altered). SELECTED_OUTPUT is used to produce a file of all the data that appear in table 22 and that were used to construct figure 6. SELECTED_OUTPUT specifies that the log of the activities of potassium ion, hydrogen ion, and silicic acid; the saturation indices for gibbsite, kaolinite, K-mica, and K-feldspar; and the total amounts in the phase assemblage and mole transfers for gibbsite, kaolinite, K-mica, and K-feldspar will be written to the file ex6A-B.sel after each calculation. The definitions for SELECTED_OUTPUT remain in effect for all simulations in the run, until a new SELECTED_OUTPUT data block is read, or until writing to the file is suspended with the identifier -selected_output in the PRINT data block.

Table 22. --Selected results for example 6A

[Simulation refers to labels in the input data set for example 6A. Negative mole transfers indicate dissolution, positive mole transfers indicate precipitation. Point on graph refers to labeled points on figure 6]



Log activity


Mole transfer, micromoles


xxxxxxx Saturation index xxxxx

Point on graph

H +


H 4 SiO 4




































































































Simulation 6A1 allows K-feldspar to react until equilibrium with gibbsite is reached. This is set up in EQUILIBRIUM_PHASES input by specifying equilibrium for gibbsite (saturation index equals 0.0) and an alternative reaction to reach equilibrium, KAlSi3O8 (the formula for K-feldspar). A large amount of K-feldspar (10.0 mol) is present to assure equilibrium with gibbsite can be obtained. Kaolinite, K-mica, and K-feldspar are allowed to precipitate if they become saturated, but they can not dissolve because they were given zero initial moles in the phase assemblage. The amount of reaction that is calculated in this simulation is the dissolution of precisely enough K-feldspar to reach equilibrium with gibbsite, possibly including precipitation of one or more of the other minerals. No gibbsite will dissolve or precipitate; the alternative reactant (KAlSi3O8) will dissolve or precipitate in its place. Simulations 6A2-6A4 perform the same calculations for kaolinite, K-mica, and K-feldspar. At other temperatures or using other minerals, it may be that a target phase is undersaturated regardless of the amount of the alternative reaction that is added because the phase is unstable relative to other phases. In this case, the numerical method will find the amount of the alternative reaction that produces the maximum saturation index.

Figure 6. --Phase diagram for the dissolution of K-feldspar (microcline) in pure water at 25 o C showing stable phase-boundary intersections (example 6A) and reaction paths across stability fields (example 6B). Diagram was constructed using thermodynamic data for gibbsite, kaolinite, K-mica (muscovite), and microcline from Robie and others (1978).

Selected results for simulations 6A1-6A4 are presented in table 22 and are plotted on figure 6 as points A, B, D, and F. The stability fields for the phases, which are based on the thermodynamic data, are outlined on the figure but are not calculated by the modeling in these simulations. At point B gibbsite starts to be transformed into kaolinite, a reaction which consumes Si. The reaction path must follow the gibbsite-kaolinite phase boundary to some intermediate point C until all gibbsite is converted, and then the path crosses the kaolinite field to point D. Similarly, there is a point E on the kaolinite-K-mica phase boundary, where the reaction path starts crossing the K-mica field to point F. Simulations 6A5 and 6A6 (table 21) solve for these two points. In simulation 6A5, point C is calculated by allowing K-feldspar to dissolve to a point where kaolinite is at saturation and is present in the phase assemblage, while gibbsite is at saturation, but not present in the phase assemblage. Likewise, simulation 6A6 solves for the point where K-mica is at saturation and present in the phase assemblage, while kaolinite is at saturation, but is not present in the phase assemblage. Assigning an initial amount of 1 mol to kaolinite in 6A5 and K-mica in 6A6 is arbitrary; the amount must only be sufficient to reach equilibrium with the mineral.

A simpler approach to determining the reaction path is to react K-feldspar incrementally, allowing the stable phase assemblage among gibbsite, kaolinite, K-mica, and K-feldspar to form at each point along the path. The only difficulty in this approach is to know the appropriate amounts of reaction to add. From points A and F in table 22, K-feldspar dissolution ranges from 0.03 to 190.9 mmol. In part 6B (table 21) a logarithmic range of reaction increments is used to define the path (solid line) across the phase diagram from its beginning at gibbsite equilibrium (point A) to equilibrium with K-feldspar (point F). However, the exact locations of points A through F will not be determined with the arbitrary set of reaction increments that are used in part B. The reaction path calculated by part 6B is plotted on the phase diagram in figure 6 with points A through F from part 6A included in the set of points.

Finally, in the kinetic approach, kinetic dissolution of K-feldspar is followed for varying amounts of time, while the phases gibbsite, kaolinite, and K-mica are allowed to precipitate and redissolve as the kinetic reaction proceeds. SOLUTION 1 is defined to have a small amount of dissolved K-feldspar (1e-13 moles). The solution then contains all elements related to phases in EQUILIBRIUM_PHASES, which, although not required for the program to run successfully, eliminates some warning messages.

During the integration of the reaction rates a simple dissolution rate law was assumed based on transition-state theory:


with k 1 = 1e-16 mol/cm 2 /s.

The KINETICS data block is used to enter specific data for the kinetic simulation. The stoichiometry of the kinetic reaction is the chemical formula of K-feldspar; by default the name of the rate is assumed to be a phase defined in PHASES data block and the formula of the phase is used as the stoichiometry of the reaction. It was assumed that the pristine soil contained 10 percent K-feldspar in the form of 0.1 mm cubes, and had g/cm 3 , so that A/V = 136/cm. The value of = 1.36e-11 mol/L/s is entered in the KINETICS data block with the identifier -parms (assuming that 1 kgw = 1 liter), and can be recalled as "PARM(1)" in the Basic rate definition in the RATES data block. It was assumed that the soil had already been weathered to some extent, and that only 90 percent of the initial K-feldspar was left [ -m0 2.16 and -m 1.94, where m0 indicates the initial mass (1 kg soil x 0.1 = 100 g / 278.3 g/mol = 0.359 mol/kg x 6 kg/L = 2.16 mol/L), and m the remaining mass (90 percent of 2.16 is 1.94 mol/L)]. The maximum amount of reaction for any time interval is restricted to 1e-6 moles ( -step_divide 1e-6). Time steps (seconds) are defined with the identifier -steps. INCREMENTAL_REACTIONS true causes the total time simulated to be the sum of all of the time steps (1.111111e8 seconds); each new time step starts at end of the previous time step.

Table 23. --Description of Basic program for K-feldspar dissolution kinetics and identification of phase transitions

Line number



Save initial amount of K-feldspar (1.94 moles)


Integrates K-feldspar dissolution rate over time interval given by TIME.


Identify greatest amount of K-feldspar present (least amount of reaction) at which gibbsite is saturated.


Identify greatest amount of K-feldspar present at which kaolinite is saturated.


Identify greatest amount of K-feldspar present at which kaolinite is saturated, but gibbsite is absent.


Identify greatest amount of K-feldspar present at which K-mica is saturated.


Identify greatest amount of K-feldspar present at which K-mica is saturated, but kaolinite is absent.


Identify greatest amount of K-feldspar present at which K-feldspar is saturated.


Save integrated reaction.


End of Basic program


Subroutine for saving values for phase transitions. If amount of K-feldspar is greater than current saved value for the index i , save amount of K-feldspar, cumulative time, log activity ratio of potassium ion divided by hydrogen ion, and log activity of silicic acid.

The rate for K-feldspar dissolution is defined in the form of Basic statements in the RATES data block. To demonstrate some of the features of the Basic interpreter, the Basic program also identifies and saves information at phase transitions, which is printed at the end of the run via USER_PRINT. The accuracy of locating a phase transition is determined by the user-definable accuracy of the integration. A small tolerance ( -tol), a large -step_divide that is greater than 1 (initial time interval will be divided by this number), or a small -step_divide that is less than 1 (specifies maximum moles of reaction) will force smaller time intervals and more accurate identification of phase transitions. In this example, -step_divide is set to 1e-6, which limits the maximum amount of reaction for any time interval to be less than 1 micromole. Thus, the amount of reaction to reach a phase transition should be identified with an accuracy of one micromole. However, limiting the amount of reaction requires smaller time intervals during the integration and, consequently, more time intervals to complete the integration, which increases the CPU time of the run.

The functions of the different parts of the Basic program are described in table 23. The functions PUT, GET, and EXISTS are used to manipulate data in static, global storage. The subscripts used in the PUT statement identify a datum uniquely. EXISTS can be used to determine if a datum with a given set of subscripts has been stored, and GET is used to retrieve data that have been stored. Once a datum has been stored with PUT, it exists for the remainder of the run, unless it is overwritten with another PUT statement with the same set of subscripts. Data stored with PUT can be retrieved by any Basic program, including those defined in RATES, USER_PRINT, and USER_PUNCH. In this example, data are stored by the RATES Basic program and the USER_PRINT Basic program retrieves the data and prints a summary of the phase transitions. While the RATES program is run many times during the kinetic integration of a time step (integration over many time intervals may be necessary for the required accuracy), the USER_PRINT program is run once at the end of each integration time step.

Table 24 gives the phase transitions encountered by the end of the last time step of example 6C. For each phase transition, the time at which the phase transition occurred, the total amount of K-feldspar that has reacted, and the coordinates of the transition on figure 6 are given. Although the values in table 24 are approximate, the amount of K-feldspar and the coordinates of the transition can be compared to table 22. As expected, the approximate mole transfers to reach the phase transitions are within 1 micromole of the values in table 22.

Table 24. --Phase transitions identified by the RATES Basic program and printed to the output file by the USER_PRINT Basic program in example 6C, which simulates the kinetic dissolution of K-feldspar

----------------------------------User print-----------------------------------
        Transition                 Time   K-feldspar      LA(K/H)   LA(H4SiO4) 
A: Gibbsite                 1.10000E+03  1.40477E-07  3.57546E-01 -6.37631E+00 
B: Gibbsite  -> Kaolinite   1.74339E+05  2.20636E-06  2.56086E+00 -5.19500E+00 
C: Gibbsite  -> Kaolinite   2.39290E+05  3.02835E-06  2.83517E+00 -5.19434E+00 
D: Kaolinite -> K-mica      1.58687E+06  2.00699E-05  4.40800E+00 -4.46585E+00 
E: Kaolinite -> K-mica      2.59719E+06  3.27914E-05  4.41031E+00 -4.25093E+00 
F: K-mica    -> K-feldspar  4.78404E+07  1.90721E-04  5.48792E+00 -3.55397E+00 

The SELECTED_OUTPUT data block specifies that a new selected-output file will be used for this simulation, ex6C.sel , and all printing to the selected-output file is eliminated ( -reset false). The USER_PUNCH data block causes two columns to be written to each line of the selected-output file, the log of the ratio of the activities of potassium ion to hydrogen ion and the log activity of silicic acid. The data will be written after each time step has been simulated ( -steps, KINETICS data block). The following table shows the results written to ex6C.sel .

Table 25. --Results written to the selected-output file by the USER_PUNCH Basic program in example 6C, which simulates the kinetic dissolution of K-feldspar

["\t " indicates a tab character]

" pH+log[K]"\t  "log[H4SiO4"\t  
-6.00016E+00\t  -1.25235E+01\t  
-1.89754E+00\t  -8.42120E+00\t  
-8.53163E-01\t  -7.37985E+00\t  
 3.57546E-01\t  -6.37631E+00\t  
 2.18230E+00\t  -5.38178E+00\t  
 4.13600E+00\t  -4.60047E+00\t  
 5.26137E+00\t  -3.71874E+00\t  
 5.48836E+00\t  -3.55363E+00\t  

| Next|| Previous || Top |