In this example, the precipitation of phases as a result of incongruent dissolution of K-feldspar (microcline) is investigated. Only the four phases originally addressed by Helgeson and others (1969)--K-feldspar, gibbsite, kaolinite, and K-mica (muscovite)--are considered. The thermodynamic data for the phases (table 21, PHASES keyword) are derived from Robie and others (1978) and are the same as for 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 (simulation 6A, table 21), the reaction path can be calculated incrementally (simulation 6B, table 21), or the reaction path can be calculated as a kinetic process (simulation 6C, table 21). 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 by 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 here. 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 obtained easily. 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, example 6 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. It is assumed that only gibbsite, kaolinite, or K-mica can form, and that these phases will precipitate reversibly if they reach saturation. That is, the phases that precipitated at the beginning of the reaction are allowed to redissolve, if necessary, to maintain equilibrium as the reaction proceeds.
The input file (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 file replaces any previous definitions for the duration of the run (the database file is not altered). SELECTED_OUTPUT is used to write a file of data that was used to produce table 22, and USER_GRAPH produces the chart shown in figure 7. SELECTED_OUTPUT specifies that the log of the activities of the 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.
A |
||||||||||||||
B |
||||||||||||||
D |
||||||||||||||
F |
||||||||||||||
C |
||||||||||||||
E |
||||||||||||||
Simulation 6A1 allows K-feldspar to react until equilibrium with gibbsite is reached. This is set up in the 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 ensure that equilibrium with gibbsite can be obtained. Kaolinite, K-mica, and K-feldspar are allowed to precipitate if they become saturated (which does not occur in this part of the simulation), but they cannot dissolve because they were given zero initial moles in the phase assemblage. The amount of reaction is the dissolution of precisely enough K-feldspar to reach equilibrium with gibbsite. 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, a target phase may remain undersaturated regardless of the amount of the alternative reaction that is added because the phase is unstable relative to the other defined phases. If this were the case, the numerical method would find the amount of the alternative reaction that produces the maximum saturation index for the target phase.
Selected results for simulations 6A1-6A4 are presented in table 22 and are plotted on figure 7 as points A, B, D, and F. The stability fields for the phases, which are based on the thermodynamic data, are calculated and plotted with USER_GRAPH as explained later. At point B, gibbsite starts to be transformed into kaolinite, a reaction that consumes Si. From this point, the reaction path follows the gibbsite-kaolinite phase boundary until all gibbsite is converted (point C), and then 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 the two points C and E. In simulation 6A5, point C is calculated by specifying that kaolinite is present at equilibrium, while K-feldspar dissolves until gibbsite is at saturation, but with zero concentration 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 with zero concentration (point E). Assigning an initial amount of 1 mol to kaolinite in 6A5 and to K-mica in 6A6 is arbitrary, but the amount must 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 to F in figure 7, K-feldspar dissolution ranges from 0.03 (simulation 6A1, table 22) to 190.9 mmol (simulation 6A4, table 22). In simulation 6B (table 21) a logarithmic range of reaction increments is used to define the path (solid line) across the phase diagram (fig. 7) 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 simulation 6B (table 21). The reaction path calculated by simulation 6B is plotted on the phase diagram in figure 7 with points A through F from simulation 6A (table 21) included in the set of points.
Finally, in the kinetic approach in simulation 6C (table 21), kinetic dissolution of K-feldspar is followed with 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 (1×10 -13 mol). 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:
, (20)where R is the rate, mol cm -2 s -1 (mole per square centimeter per second); k 1 is the rate constant, equal to 1×10 -16 mol cm -2 s -1 ; m is the amount of kinetic reactant, mol; m 0 is the initial amount of kinetic reactant, mol; IAP is the ion activity product; and K is the equilibrium constant.
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 the 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.36×10
-11
mol L
-1
s
-1
(mole per liter per second) 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 × 0.1 = 100 g / 278.3 g/mol = 0.359 mol/kg × 6 kg/L = 2.16 mol/L), and
m
is 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 1×10
-6
mol (
-step_divide
1e-6). Time steps (s) are defined with the identifier
-steps
. INCREMENTAL_REACTIONS
true
causes each new time step to start at the end of the previous time step, so that the total time is the sum of all the time steps (= 1.111111×10
8
s).
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 by using 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 simulation (6C), -step_divide is set to 1×10 -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 will be identified with an accuracy of 1 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 purpose of each part of the Basic program is described in table 23. The functions PUT, GET, and EXISTS are used to manipulate data in static, global storage. The arguments used in the PUT function identify a number and the storage location. For example, PUT(M,i) will place the value of “M” in the location identified by the variable i. EXISTS can be used to determine if a storage location contains a number, and GET is used to retrieve data that have been stored. For example, IF EXISTS(i) then m2 = GET(i) will assign the number in location i to m2 if i contains a number, and do nothing if it does not. Once a number has been stored with PUT, it exists for the remainder of the run, unless it is overwritten with another PUT statement in the same location (that is, with the same set of subscripts). Data stored with PUT can be retrieved by any Basic program defined in RATES, USER_GRAPH, USER_PRINT, and USER_PUNCH. In this simulation (6C), 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. Whereas the RATES program is run many times during the kinetic integration of a time step, the USER_PRINT program is run only once at the end of each integration time step.
See Phase transitions identified by the RATES Basic program and printed to the output file by the USER_PRINT Basic program in simulation 6C, which simulates the kinetic dissolution of K-feldspar. gives the phase transitions encountered by the end of the last time step of simulation 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 7 are given. Although the values in table 24 are. approximate, the amount of K-feldspar and the coordinates of the transition can be compared with table 22 As expected, the approximate mole transfers to reach the phase transitions in table 24 (column 3) are within 1 micromole of the values in table 22 (column 2).
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 are written after each time step has been simulated ( -steps , KINETICS data block). See Results written to the selected-output file by the USER_PUNCH Basic program in simulation 6C, which simulates the kinetic dissolution of K-feldspar. shows the results written to ex6C.sel .
The phase boundaries in figure 7 are plotted with USER_GRAPH at the end of the file. First, a dummy phase “K_H” is defined that will provide the activity ratio of K+ and H+ in the solution, which is plotted on the y-axis. The -log_k of zero for the reaction is the default value in PHREEQC and may be omitted. The saturation index of a solution for this phase is given by log([K+] / [H+]), and can be obtained with SI("K_H") in a Basic program. Next, four solutions are defined for the points that delimit the phase boundaries, going from top to bottom and from left to right. USER_GRAPH connects the points with a black line, without symbols, and without entry in the legend (no - headings). After the END, another four solutions are defined for the phase boundaries from top to bottom and from right to left, and another curve is drawn by redefining line 10 in USER_GRAPH (in this way avoiding connecting the last bottom-right point with the new top-right point). Note that the text is compacted in this part of the input file by concatenating lines with semicolons, “;”.
