Environ. Sci. Technol. 2004, 38, 3641-3648
Predominance and Mineral Stability Diagrams Revisited D A V I D G . K I N N I B U R G H * ,† A N D DAVID M. COOPER‡ British Geological Survey and Centre for Ecology and Hydrology, Wallingford, Oxon OX10 8BB, U.K.
Predominance and mineral stability diagrams are most useful if they are based on a full chemical speciation. The diagrams can then include a wide range of reactions of importance in environmental chemistry including adsorption and ion exchange reactions. They can also take into account the possible interactions between a large number of components and are able to model changes at a fixed total concentration of each component rather than under the commonly used but less realistic assumption of a fixed activity. A “hunt-and-track” algorithm is described that enables predominance and mineral stability diagrams to be constructed using any general-purpose speciation program. Examples are given using PHREEQC and Orchestra.
Introduction Predominance diagrams show the dominant or major species in solution as a function of two or more master variables, often the activities or activity ratios of the major controlling species present. A mineral stability diagram is somewhat similar but shows the boundaries between the stability regions of different minerals (i.e., where the mineral begins to precipitate). Garrels and Christ (1) popularized the use of such diagrams in geochemistry following the extensive earlier work of Pourbaix in corrosion science (2). Pourbaix (2), and later Garrels and Christ (1), called the combination of a predominance diagram and a mineral stability diagram, a “combination” diagram. Such a diagram shows the dominant species (gas, solution, and mineral phases) that exist under a range of conditions. These conditions are specified by the systematic variation of two master variables represented by the x- and y-axes. The Eh-pH or pe-pH diagram is perhaps the best known such diagram. Such diagrams now feature in virtually every textbook on geochemistry and aquatic chemistry (3-6). The basic principles involved in their calculation are described in many textbooks. One of the clearest descriptions is in Morel and Hering (7). Many software packages are available for calculating such “classical” predominance and stability diagrams, but probably the most elegant and widely used in the field of geochemistry is the Act2 package in The Geochemist’s Workbench (8). Given a range of activities for the x and y master variables plus a specified component, such as iron, and the activities of other important ions, Act2 is able to determine and plot the dominant species and mineral phases that are stable over the specified region. Act2 is interactive and easy to use. Unfortunately, many other excellent speciation programs do not currently include the * Corresponding author phone: +44(0) 1491 69 2293; fax: +44(0) 1491 69 2345; e-mail:
[email protected]. † British Geological Survey. ‡ Centre for Ecology and Hydrology. 10.1021/es034927l CCC: $27.50 Published on Web 05/22/2004
2004 American Chemical Society
ability to calculate such diagrams. This includes PHREEQC (9), MINTEQA2, MINEQL, and Orchestra (10). It is commonly necessary to invoke certain assumptions when constructing classical predominance diagrams, namely, that (i) the diagrams are calculated on the basis of a fixed activity of the main component and known activities of all minor species; (ii) the activity of all solids is equal to unity; and (iii) there is usually a limit to the number of reactive components considered. The resulting boundaries are often linear when plotted on a log-log diagram, with a slope of 1:1, 1:2, etc. depending on the stoichiometry of the dominant reaction defining the boundary. Complications arise when the form of the main species varies in the chosen x-y space, but this can be allowed for by creating a series of subdiagrams and assembling the main diagram from these, as in Act2 (8). These assumptions may not be too serious given that the diagrams are usually only used in a semiquantitative way. However, what is probably of greater concern, and what often limits their practical usefulness, is that they omit certain types of reaction completely, namely, reactions involving adsorption, ion exchange, and solid solution formation. They also omit complexation reactions that fall outside the normal computational scheme such as those describing metalhumic interactions. These reactions can dominate the behavior of trace metals in environmental systems and should therefore not be ignored. Because of the rather unspecific nature of many adsorption and ion exchange reactions and the importance of competitive interactions (the essence of ion exchange), these reactions are wide-ranging and common. All possible interactions must be considered because of their possible impact on the final speciation. This requires input of the complete water chemistry and a full speciation. It would be difficult to include the additional features identified above in the classical scheme because the simplifications invoked during their construction no longer apply. For example, in surface complexation models that consider electrostatic interactions between the adsorbed ions, the Boltzmann factor appears as an activity coefficient. This is often far from unity and varies strongly with environmental factors such as pH, ionic strength, etc. It is a key part of the surface complexation calculation and must be treated as an unknown to be calculated by iteration.
New Approaches Based on a Full Chemical Speciation Rather than adopting the classical “analytical” approach, a numerical approach can be adopted in which no attempt is made to produce a simple set of equations defining the boundaries. Rather a full speciation is carried out using any general purpose speciation program, and then the dominant species is simply calculated from the results. Such an approach is far more flexible and removes many of the limitations inherent in the classical approach. This approach has been adopted by programs such as MEDUSA, FactSage, and CSP (11-14). However, these programs and their thermodynamic databases tend to reflect their origins in industrial chemistry rather than in environmental chemistry. They do not include adsorption reactions or many of the other reactions of interest to environmental chemists. Their thermodynamic databases are also often proprietary and hidden from the user. Therefore, we have explored the use of more environmentally orientated programs for calculating predominance diagrams. The minimum that is required is that the speciation program accepts a given a set of solution conditions and constraints and returns the concentrations of the dominant species present (gas, solution, adsorbed, and mineral). VOL. 38, NO. 13, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3641
FIGURE 1. Hunt-and-track strategy for locating the boundary between two fields. (a) After locating a change in dominance on a domain boundary, cells are explored in the order a, b, c, ... . The numbers indicate the order and location where speciation calculations are undertaken, and the color of the symbols at the grid intersections indicates the dominant species returned. (b) Method of linear interpolation used to establish the location of the boundary as it enters cell c (see text and eq 1). Communication can be either through files or if the speciation program permits by tighter integration through a more direct interface. We have adopted the former approach using the PHREEQC (9) and Orchestra (10) programs exactly as distributed. We also tested an early version of iPHREEQC, a special version of PHREEQC that is callable from other programs. PHREEQC is a well-supported and well-documented public domain speciation and transport program that is widely used by geochemists. It includes options for solution complexation, gas solubility, mineral precipitation/dissolution, solid solution formation, adsorption, ion exchange and surface complexation, and reaction kinetics and includes an optional GUI program to help construct the input files and to view the output files. Its associated thermodynamic databases are stored in text format and can be readily inspected, edited, and extended by the user. PHREEQC uses either the Davies equation or an extended form of the DebyeHu ¨ ckel equation for ionic strength corrections for charged species. The Dzombak and Morel (15) diffuse layer model is included for surface complexation modeling. PHREEQC runs as a batch program using separate input, output, and database files. We wrote a small driver program to run PHREEQC in “hunt-and-track” mode. This repeatedly writes the input file in response to requests from the tracking program and reads the results from a PHREEQC selected output file. Orchestra (10) is a novel and flexible Java program for creating sophisticated geochemical models including a wide range of adsorption and surface complexation models. It too can be run in batch mode via input and output files. Hunt-and-Track Algorithm. Our strategy for locating the domain boundaries makes use of the fact that all domain and axis boundaries are interconnected (i.e., there are no “islands”). We assume that, presented with values for the two master variables, the speciation program can return the dominant species. In principle, it is a simple matter to build up a complete diagram by calculating the dominant species on a grid. If the grid is rectangular, the results can be rendered directly, pixel by pixel. This is the approach recently adopted by Solgaswater for Windows (16). However, it is more elegant to locate and parametrize the field boundaries. This can be 3642
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 13, 2004
done using raster-to-vector, GIS, or custom software. This is the approach adopted by Beverskog and Puigdomenech (17) and implemented by the MEDUSA (14) program. This “brute force” approach is rather inefficient since many of the points evaluated will be far from a boundary and will therefore contribute little useful information. It is also difficult to guarantee locating all of the fields without a fairly large computational effort. Therefore we have adopted a different strategy which we term the hunt-and-track approach. This concentrates on locating and tracking the various boundaries and puts no effort into computations far from the boundaries. It requires no initial grid of points. The algorithm is apparently similar to that sometimes used to locate “zero phase fraction lines” in phase diagrams. Our hunt-and-track algorithm is based on a fixed regular grid over the domain of interest defined by the x- and y-axes. The domain is usually rectangular, with its boundaries (axes) also being grid lines. The algorithm assumes that at least one predominance boundary actually crosses an axis (if none does, there must be just a single dominant species). The first task is to identify such a crossing point by searching along the domain boundaries or axes. Starting at a selected corner of the domain, the dominant species is identified. The next point for computation is the nearest grid point along the axis, moving in a selected direction. If this has the same dominant species, the algorithm moves on to successive grid points on the axis until a change is found. This is the “hunting” part of the algorithm. Once a point of change has been identified, the algorithm moves into “tracking” mode, following the predominance boundary within the domain. Figure 1 shows an example with a predominance boundary (shown as a dashed line), whose precise location is unknown, crossing the y-axis of a rectangular domain. Grid points where speciation calculations are undertaken are identified in Figure 1a by numbered circles. The numbering indicates the order in which the grid points are visited. In this example, the first grid point in the sequence is on the y-axis where the algorithm is in hunting mode. It continues through points 2 and 3, where there is no change in dominant species, reaching point 4 where a change is identified. There must therefore be an intersection of the
y-axis with the predominance boundary between points 3 and 4. The tracking mode now begins with the dominant species calculated at points 5 and 6, which complete rectangle a on the grid. The predominance boundary must exit through one of the remaining three sides of this rectangle. This is immediately identified as the side linking points 5 and 6. This exit side for rectangle a becomes the entry side for rectangle b, the second to be constructed. Calculation of the dominant species at points 7 and 8 shows that the exit side for rectangle b is the side linking points 5 and 7. Rectangle c is built on this side, with the exit side linked by points 7 and 10, on which rectangle d is constructed. This technique tracks a sequence of grid squares through which the predominance boundary runs. It does not provide coordinates through which the line passes. These can be approximated by linear interpolation along each exit or entry side. Four values are required to carry out this interpolation (Figure 1b). These are the concentrations of the dominant and subdominant species at the end points of the side. Denote the end points (x1, y1) and (x2, y2) and the dominant and subdominant log concentrations as d1, s1, d2, and s2. Then the linearly interpolated approximate location (xb, yb) of the predominance boundary as it crosses the line joining (x1, y1) and (x2, y2) is
xb ) x1 +
d 1 - s1 (x - x1) d2 - s2 + d1 - s1 2
d1 - s 1 yb ) y1 + (y - y1) d2 - s2 + d1 - s1 2
(1)
On a rectangular grid, either x1 ) x2 or y1 ) y2, so one or other of the equations will always be degenerate, with either xb ) x1 or yb ) y1. It is straightforward to extend this concept to three components. This improves the location of triple points. Such an interpolation approach no longer applies to boundaries involving minerals. Here the boundary is set to halfway between the end points. The successive construction of exit/entry lines and rectangles on the fixed grid continues until either it exits from the domain or a junction on the predominance boundary is identified. A junction is indicated when more than two species are identified as dominant among the four corners of a rectangle. If, for example, there are three dominant species, then these define two exit sides from a single rectangle. The two boundary lines from this rectangle are then tracked in turn, with appropriate flagging of the necessary sequence of operations. For more complex examples, a predominance boundary may return to a previously identified junction. By suitable housekeeping, the algorithm identifies when this occurs. Finally, the line segments are assembled into polygons. The domain boundary should be exhaustively searched in hunting mode to ensure that no predominance boundaries are missed. Stepping sequentially through every domain boundary grid point is recognized as inefficient, and time can be saved by adopting a binary search. However, the subdominant species must also be monitored to ensure that no boundary is missed. If the grid is too coarse, the location of junctions may be poorly estimated. For a fixed grid algorithm of this kind there will be a tradeoff between computational efficiency and accuracy. A grid spacing equivalent to a resolution of 1/500th of the domain boundary normally gives smooth and reliable curves. The advantage of using a fixed grid is that quite complex curves may be easily tracked. However, the scheme does not take full account of the known characteristics of predominance boundaries inherent in the underlying chemical model, and this can lead to some inefficiency. A similar tracking procedure can
be adopted when the stability criterion (see below) is used or when an artificial boundary is added as a result of a userdefined constraint (e.g., PH2 e 1 atm). Definitions of Predominance and Mineral Stability. In classical mineral stability diagrams, the position of the mineral-solution boundary depends on the assumed activity of the solution components at the boundary. Garrels and Christ, following Pourbaix, defined the boundary as the point where the “sum of the activities of the ions in equilibrium with the solid exceeds some chosen value”. They chose 10-6 as a default value on the basis that if it is less than this value, the solid will tend to behave as an immobile constituent in the environment. In the full speciation approach, the activity at the boundary is determined by the system specification and the speciation calculation. It will normally be specified by either a fixed total amount of the various components or by a fixed activity or fugacity (for gases). There is no need for excessive simplification, and so the diagrams can be drawn for systems of particular interest. The only decision is whether to draw the predominance boundary or the stability boundary (i.e., at mineral-solution boundaries, it is either the point where the dominant species changes (ci > cj) or the point where mineral precipitation begins (ci > 0)). It is useful to view predominance in terms of the natural chemical hierarchy that exists within chemical systems (10). This can be viewed from different levelssthe “phase” level (gas, solution, adsorbed, solid) and the species level (Fe3+, FeOH2+, ...). Normally, the diagrams are drawn at the species level. The definition of “species” becomes more ambiguous when an adsorbed phase is considered since a common definition of an adsorbed species is in terms of the type of binding site. However, this would make the diagrams very sensitive to the adsorption model adopted, and so it is probably more useful to treat all adsorbed entities for a given adsorbent as a single species for the purposes of ranking and display. If the adsorbed speciation is of particular interest, then treating each entity separately can show this.
Applying the New Approach The following examples are intended to illustrate some of the important features observed during the calculation of the predominance and stability diagrams using the new approach. It is not the intention here to assess critically the databases used, the real-world relevance of the reactions included, or the activity coefficient models used. Iron Hydrolysis and Reduction (Fe-H2O-Cl). The following examples (i) demonstrate how even a “simple” system can give rise to complex diagrams when viewed using the full chemical speciation approach and (ii) illustrate how the various choices (e.g., the predominance vs stability option, database used) can influence the final output. The hydrolysis and reduction of Fe3+ as a function of pH and redox status is an important and complex reaction that is of considerable environmental significance. It provides a good test of the hunt-and-track algorithm. The following results were based on the llnl.dat database (revision 1.12) distributed with PHREEQC. This is the most comprehensive of the four standard PHREEQC databases and should give speciation results essentially equivalent to those from EQ3/6 and The Geochemist’s Workbench. Most of the thermodynamic data for the hydrolyzed species are taken from the compilation of Baes and Mesmer (18). The extended DebyeHu ¨ ckel equation with A, B, and B-dot parameters is used for estimating the activity coefficients. This also includes an ionspecific ion size parameter, a. Other choices of thermodynamic data are possible (17). The pH of a solution containing 0.01 mol/kgw FeCl3 and 0.1 mol/kgw NaCl (kgw ) kg of water) initially at pH 1.95 was adjusted over the pH range 2-12 by titrating with solid NaOH. The redox status was adjusted by equilibrating with various VOL. 38, NO. 13, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3643
FIGURE 2. Predominance (a, b) and stability (c, d) diagrams for the Fe-H2O-Cl system calculated with PHREEQC using the llnl.dat thermodynamic database but excluding the crystalline iron(III) oxides and magnetite (FeT ) 0.01 mol/kgw, 0.1 mol/kgw NaCl, 25 °C). partial pressures of O2(g) (log PO2(g) -85 to 0) (partial pressures can be equated with fugacities for an ideal gas). Constraints of PO2 e 0.21 atm and PH2 e 1 atm were included based on the limits imposed by the decomposition of water. Charge balance was maintained by adding Cl- as necessary. Initially only a small subset of minerals (Fe(OH)3(a), Fe(OH)2, and FeO) were allowed to precipitate. The more stable (and insoluble) minerals such as hematite, goethite, and magnetite were excluded. Calculations were made at 25 °C using both predominance and stability options. Hydrous ferric oxide (Fe(OH)3(a) or HFO) was found to be stable over much of the region with FeO being stable under strongly reducing conditions and pH > 8 (Figure 2). The predominance and stability fields were broadly similar with, as expected, slightly larger fields for the two solid phases under the stability option. Close-ups of the pH 3.0-3.5 region show more subtle differences: the predominance diagram shows that the FeOH2+ field is split into two by the Fe3(OH)45+ field. There is also a small triangular FeOH2+ field (at pH 3.32, log PO2 -22.4) to the right of the main FeOH2+ field. This field is not seen in the corresponding stability diagram as it is covered by the enlarged Fe(OH)3(a) field. Note also the strongly curved boundary between the Fe2(OH)24+ field and the main FeOH2+ field. Such curved boundaries are not seen in classical stability diagrams. The same calculations were repeated but magnetite (Fe3O4) was included as a possible mineral phase. Magnetite has a broad range of stability under reducing and nonacidic conditions and replaces Fe(OH)3(aq) and FeO as the most stable mineral under these conditions (Figure 3). In the pH range 2.3-2.8, the predominance diagram shows a very small FeOH2+ field at pH 2.6, log PO2 ) -22. This extends to higher 3644
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 13, 2004
pH and PO2 as a narrow, parallel field between the Fe2(OH)24+ and magnetite fields. These small wedge-shaped fields disappear when the stability criterion is used. A pH section at log PO2 ) -13 shows how the concentrations of the various Fe species vary between pH 2.9 and 3.0 (Figure 4). Magnetite precipitation starts at about pH 2.92 and becomes the dominant form of Fe at about pH 2.95. The vertical arrows denote the locations of the boundary defined according to the predominance and stability options. If the highly stable mineral hematite (R-Fe2O3) is allowed to precipitate, then this dominates much of the stability diagram (Figure 5a). The diagrams are clearly highly dependent on the thermodynamic data used and the choice of minerals allowed to precipitate. The alternative wateq4f.dat database (revision 2.5), also distributed with PHREEQC, includes the mineral Fe(OH)2.7Cl0.3. If this database is used, this mineral becomes dominant over much of the domain (Figure 5b). Carbonate Equilibria (Ca-Mg-Zn-CO3). Dissolved inorganic carbon species are an important component of all natural waters. In open systems, their concentrations are controlled by the pH, PCO2, and the presence of metal carbonate complexes and mineral phases. At high pH, the carbonate concentration in equilibrium with the atmosphere is very high, and so it is useful to be able to remove such unrealistic regions from the diagram. This example illustrates this. Predominance diagrams have been calculated for the system Ca-Mg-Zn-CO3 based on an initial solution containing 5 mmol/kgw of calcium and zinc chlorides and 2 mmol/kgw MgCl2 at pH 7 and subject to varying PCO2 (by the addition of CO2(g)) and pH (by the addition or subtraction of NaOH). The area where the total carbon concentration
FIGURE 3. Same as for Figure 2 but allowing magnetite to precipitate (FeT ) 0.01 mol/kgw, 0.1 mol/kgw NaCl, 25 °C). system oxidizing. The pH-log PCO2 predominance diagrams were calculated for the various components. As expected, the calcium diagram (Figure 6) shows that Ca2+ is dominant below pH 7-8 after which calcite becomes the major form of Ca. About one-third of the diagram is blank because this is where CT > 0.3 mol/kgw. Dolomite is stable under broadly similar conditions to calcite. Zn2+ tends to form the minerals zincite (ZnO), hydrozincite, and Smithsonite (ZnCO3) above pH 6-7 depending on the PCO2. At high pH and high PCO2 values, the aqueous species Zn(CO3)22- can dominate the Zn speciation. The carbonate diagram shows the overlapping contributions from all of these reactions. At low PCO2 values and high pH, there is an excess of carbonate in the system and CO32- and NaCO3- predominate. The Na is derived from the NaOH used to raise the pH. As the PCO2 is increased, the total carbonate constraint comes into effect at an increasingly low pH as expected.
FIGURE 4. Variation in iron concentration and the dominant iron species near the narrow field boundary at log PO2 ) -13 seen in Figure 3. The name of the dominant species is given at the bottom of the diagram. The FeOH2+ field completely disappears under the stability criterion. exceeds 0.3 mol/kgw has been blanked out (the ion activity model used is likely to fail here). The wateq4f.dat database was used, but data for the formation of hydrozincite, Zn5(OH)6(CO3)2, was added (19):
Zn5(OH)6(CO3)2 + 10H+ ) 5Zn2+ + 2CO2 + 8H2O log K ) 45.0; ∆H ) -256.5 kJ No minerals were suppressed and PO2 was maintained at atmospheric concentrations (log PO2 ) -0.678) to keep the
Acid Mine Water. The following example shows how the new approach enables pe-pH diagrams to be calculated for complex systems. It includes adsorption linked to two mineral phases and a CO2 constraint, albeit imposed slightly differently from the example above. The diagram was constructed for a sample of Fe-, As- and SO4-rich mine water (WJ1) from Wheal Jane Mine, Cornwall, U.K. (20). The pH was 3.5, and the estimated pe was 1.9. The initial composition was (in mg/kgw) as follows: Cl 179, F 44, SO4 1390, Ca 191, Mg 43, Na 93, K 12, Al 25, Si 11.0, Sr 1.87, Ba 0.052, Li 2.7, Fe 346, Mn 19.7, As 2.1, and Zn 125. These concentrations were fixed, and the pH was adjusted by the addition (or subtraction) of NaOH. Na was supplied by equilibration with a fictitious halite (NaCl) phase (SI ) -6.34) to ensure that sufficient Na was present to achieve the lowest pH values. The redox status was controlled by titration with O2 gas. The CO2(g) partial pressure was maintained at its atmospheric concentration (log PCO2 ) -3.5) subject to the constraint that the total amount of carbon that could be supplied should not exceed VOL. 38, NO. 13, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3645
FIGURE 5. Stability diagrams for the Fe-H2O-Cl system including (a) magnetite and hematite (llnl.dat database) and (b) magnetite, hematite, and Fe(OH)2.7Cl0.3 (wateq4f.dat database) as possible stable minerals (FeT ) 0.01 mol/kgw, 0.1 mol/kgw NaCl, 25 °C).
FIGURE 6. Ca, Mg, Zn, and C predominance diagrams for the Ca-Mg-Zn-CO3-Na-Cl system as a function of pH and log PCO2. Areas where the total carbon concentration in the system exceeds 0.3 mol/kgw have been omitted. 0.01 mol/kgw. Where the amount of CO2(g) required exceeded this amount, only 0.01 mol/kgw was added. This meant that the final PCO2 was lower than the target value. The wateq4f.dat database was used, which included the revised set of arsenic data of Nordstrom and Archer (21). Hydrozincite was also included. Fe(OH)3(a) () hydrous ferric oxide or HFO) was the only iron(III) oxide mineral allowed to precipitate. Many other kinetically unfavorable minerals, including most of the silicates, were also excluded. When present, the HFO adsorbed cations and anions according to the Dzombak and Morel (15) diffuse layer model with default database and parameters. This takes into account the direct (competitive) and indirect (electrostatic) interactions between all the 3646
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 13, 2004
adsorbed ions. Similarly adsorption was also linked to the precipitation of Al(OH)3(a) (hydrous aluminum oxide, HAO). In the absence of a critically compiled sorption database for HAO, the default parameters for HFO were used except that the pzc was increased by 1 log unit and the specific surface area was assumed to be 100 m2/g (7800 m2/mol) rather than 600 m2/g. The predominance diagrams for Al, As, C, Ca, F, Fe, Mn, S, and Zn (Figure 7) provide varying views of the reactions occurring. Fe(OH)3(a) precipitates over much of the domain. This adsorbs many anions and cations including As(V) and As(III). Arsenic adsorbed to HFO becomes the dominant As species whenever HFO is stable (i.e., except under strongly
FIGURE 7. pe (Eh)-pH diagrams for Al, As, C, Ca, F, Fe, Mn, S, and Zn using acid mine drainage from the Wheal Jane Mine, England, as a starting point (initial pe and pH marked by the red cross) and in equilibrium with log PCO2 ) -3.5 subject to the constraint that the maximum amount of CO2 added is 0.01 mol/kgw. acidic or strongly reducing conditions). This would not be shown in classical pe-pH diagrams. Al(OH)3(a) also precipitates over a wide range of pH conditions and adsorbs As and other constituents. This is especially important for As under near-neutral and strongly reducing conditions where HAO is stable but HFO is not. Adsorbed species are also important for Zn but less so because the formation of hydrozincite and amorphous ZnO are favored at high pH and because of the relatively weak adsorption of Zn at pH < 6.5. Carbonate minerals are important for Mn, Ca, Mg, Zn, and Fe especially at high pH where carbonate concentrations are relatively high. Under highly reducing conditions, carbonate is reduced to methane and the PCH4 < 1 atm constraint forms the lower boundary to the diagram. CH4(aq) becomes the dominant C species close to this boundary under neutral to acid conditions. Manganese oxides are important under oxidizing and high pH conditions while the carbonate mineral (rhodochrosite) dominates the Mn speciation under reducing and alkaline conditions. Calcite and dolomite precipitate at high pH. Sulfate is reduced to sulfide under strongly reducing conditions and leads to pyrite (FeS2), sphalerite (ZnS), and orpiment (As2S3) formation. At nearneutral pH and strongly reducing conditions, there is a narrow field where the thioarsenite species AsS(OH)(HS)- dominates the As speciation. Fluorite (CaF2) is stable in the region pH
6.5-9 whereas soluble F, Al-F, and Fe-F species dominate the F speciation elsewhere. Cadmium Binding by Humic Acid. This example is for a class of reaction ignored in classical predominance diagrams. The binding of H+ and Cd2+ to humic acid (HA) involves strong site heterogeneity, electrostatic interactions, and a variable reaction stoichiometry. This results in strongly curved field boundaries. The NICA-Donnan model (22) is able to calculate the specific and nonspecific binding of metal ions by HA as a function of the concentration of H+, free metal ions, and HA. This model is not implemented in PHREEQC but has been implemented in Orchestra (10), and so this program was used to calculate the speciation for a total Cd concentration of 1 mmol/kgw in a 0.1 mol/kgw KNO3 background electrolyte. Thermodynamic data for the inorganic species involved have been taken from the wateq4f.dat database while proton and metal ion binding constants for HA are from Milne et al. (23, 24). Inorganic Cd species considered were Cd2+, CdOH+, Cd(OH)2, Cd(OH)3-, Cd2(OH)3+, Cd(OH)42-, and CdNO3+ in solution. Minerals considered were Cd(OH)2(s) and monteponite (CdO). At low concentrations of HA relative to the total metal ion concentration, Cd2+ dominates but as the concentration of HA is increased, Cd bound to HA becomes dominant. Above pH 8.7, the mineral species VOL. 38, NO. 13, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3647
Parkhurst, Scott Charlton, Hans Meeussen, Tony Appelo, Bob Simons, and Eric Fookes. We also thank Kirk Nordstrom for an early sight of his revised thermodynamic data for arsenic, Scott Charlton and David Parkhurst for an early version of their callable interface to PHREEQC and for other discussions, and Hans Meeussen for help with Orchestra. This paper is published with the permission of the Executive Director of the British Geological Survey (NERC).
Literature Cited
FIGURE 8. Predominance diagram for Cd in the system Cd-humic acid (HA) as a function of pH and HA concentration for a fixed total Cd concentration of 1 mmol/kgw. Cd(OH)2(s) becomes important especially at HA concentrations below 100 mg/L (Figure 8). The HA concentrations considered are low as compared with the concentrations in many soil systems where the HA concentration is typically of the order of hundreds of grams per liter.
Discussion The numerical approach to the calculation of predominance diagrams outlined here has similar advantages and disadvantages to that found when comparing numerical and analytical approaches for other applications (e.g., for contaminant transport models). Advantages are generality, flexibility, and ease of maintenance; disadvantages are computational speed, elegance, and possible numerical problems. The computational effort required is similar to that involved in a reaction path or (1D) transport calculation. We believe that the advantages of the numerical approach are overwhelming. Inclusion of all the possible processes included in the speciation model is logical and consistent. All points within the diagram can be reached by a plausible reaction path. Fixed total concentrations are more plausible than fixed activities. The numerical approach is robust and is able to model complex systems close to those of environmental interest including many of the complex interactions involved. The diagrams are also useful for comparing thermodynamic databases. At present, one of the greatest limitations is in knowing which minerals will precipitate or dissolve within the time scale being considered. As speciation software and its associated databases become increasingly sophisticated and more faithfully reflect “the real world”, so will the predominance diagrams calculated from them. Furthermore, as the complexity and range of processes considered by speciation programs increases, so does the need for simple graphical tools to reflect and highlight these processes. Predominance diagrams provide an efficient way of revealing the essential features of these complex chemical systems in an appealing way. The widespread availability of color graphics, both display and printed, further enhances their usefulness.
Acknowledgments We appreciate the dedicated efforts of the software writers who have made these calculations possible, particularly David 3648
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 13, 2004
(1) Garrels, R. M.; Christ, C. L. Solutions, Minerals and Equilibria; Harper and Row: New York, 1965. (2) Pourbaix, M. Atlas of Electrochemical Equilibria; Pergamon Press: Oxford, 1966. (3) Brookins, D. G. Eh-pH Diagrams for Geochemistry; SpringerVerlag: Berlin, 1988. (4) Fletcher, P. Chemical Thermodynamics for Earth Scientists; Longman Scientific and Technical: Harlow, 1993. (5) Stumm, W.; Morgan, J. J. Aquatic Chemistry; 2nd ed.; John Wiley: New York, 1981. (6) Pytkowicz, R. M. Equilibria, Nonequilibria and Natural Waters. Vol. II; John Wiley: New York, 1983. (7) Morel, F. M. M.; Hering, J. G. Principles and Applications of Aquatic Chemistry; John Wiley & Sons: New York, 1993. (8) Bethke, C. M. The Geochemist’s Workbench. Release 4.0. A User’s Guide to Rxn, Act2, Tact, React and Gtplot; University of Illinois: 2002. (9) Parkhurst, D. L.; Appelo, C. A. J. User’s Guide to PHREEQC (Version 2)sA Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations; Water-Resources Investigations Report 99-4259; U.S. Geological Survey: Reston, VA, 1999; wwwbrr.cr.usgs.gov/ projects/GWC_coupled. (10) Meeussen, J. C. L. Environ. Sci. Technol. 2003, 37, 1175-1182; www.meeussen.nl/orchestra/. (11) Bale, C. W.; Chartrand, P.; Degterov, S. A.; Eriksson, G.; Hack, K.; Mahfoud, R. B.; Melanc¸ on, J.; Pelton, A. D.; Petersen, S. Calphad 2002, 26, 189-228. (12) Anderko, A.; Sanders, S. J.; Young, R. D. Corrosion 1997, 53, 43-53. (13) Anderko, A.; Shuler, P. J. Comput. Geosci. 1997, 23, 647-658. (14) Puigdomenech, I. MEDUSA: Make Equilibrium Diagrams Using Sophisticated Algorithms; 2002; www.kemi.kth.se/medusa. (15) Dzombak, D. A.; Morel, F. M. M. Surface Complexation ModellingsHydrous Ferric Oxide; John Wiley & Sons: New York, 1990. (16) Hedlund, T. WinSGWsSolgaswater for Windows; Inorganic Chemistry Department, Umeå University: Umeå, Sweden, 2003; http://www.chem.umu.se/dep/inorgchem/samarbeta/ WinSGW_eng.stm. (17) Beverskog, B.; Puigdomenech, I. Corros. Sci. 1996, 38, 21212135. (18) Baes, C. F.; Mesmer, R. E. The Hydrolysis of Cations; Krieger: Malabar, 1986. (19) Preis, W.; Gamsja¨ger, H. J. Chem. Thermodyn. 2001, 33, 803819. (20) Edmunds, W. M.; Andrews, J. N.; Bromley, A. V.; Kay, R. L. F.; Milodowski, A.; Savage, D.; Thomas, L. J. Granite-water interactions in relation to Hot Dry Rock geothermal development. Investigation of the Geothermal Potential of the UK; British Geological Survey: Wallingford, 1988. (21) Nordstrom, D. K.; Archer, D. G. In Arsenic in Ground Water; Welch, A. H., Stollenwerk, K. G., Eds.; Kluwer: Boston, 2003; pp 1-25. (22) Kinniburgh, D. G.; van Riemsdijk, W. H.; Koopal, L. K.; Borkovec, M.; Benedetti, M. F.; Avena, M. J. Colloids Surf. A 1999, 151, 147-166. (23) Milne, C. J.; Kinniburgh, D. G.; Tipping, E. Environ. Sci. Technol. 2001, 35, 2049-2059. (24) Milne, C. J.; Kinniburgh, D. G.; van Riemsdijk, W. H.; Tipping, E. Environ. Sci. Technol. 2003, 37, 958-971.
Received for review August 22, 2003. Revised manuscript received January 5, 2004. Accepted March 12, 2004. ES034927L
