Uniqueness

A practical question that arises in quantitative modeling is whether the results of a modeling study are unique. In other words, is it possible to arrive at results that differ, at least slightly, from the original ones but nonetheless satisfy the governing equations and honor the input constraints? In the broadest sense, of course, no model is unique (see, for example, Oreskes et al., 1994). A geochemical modeler could conceptualize the problem differently, choose a different compilation of thermodynamic data, include more or fewer species and minerals in the calculation, or employ a different method of estimating activity coefficients. The modeler might allow a mineral to form at equilibrium with the fluid or require it to precipitate according to any of a number of published kinetic rate laws and rate constants, and so on. Since a model is a simplified version of reality that is useful as a tool (Chapter 2), it follows that there is no“correct” model, only a model that is most useful for a given purpose. A more precise question (Bethke, 1992) is the subject of this chapter: in geochemical modeling is there but a single root to the set of governing equations that honors a given set of input constraints? We might call such a property mathematical uniqueness, to differentiate it from the broader aspects of uniqueness. The property of mathematical uniqueness is important because once the software has discovered a root to a problem, the modeler may abandon any search for further solutions. There is no concern that the choice of a starting point for iteration has affected the answer. In the absence of a demonstration of uniqueness, on the other hand, the modeler cannot be completely certain that another solution, perhaps a more realistic or useful one, remains undiscovered. Geochemists, following early theoretical work in other fields, have long considered the multicomponent equilibrium problem (as defined in Chapter 3) to be mathematically unique. In fact, however, this assumption is not correct. Although relatively uncommon, there are examples of geochemical models in which more than one root of the governing equations satisfy the modeling constraints equally well. In this chapter, we consider the question of uniqueness and pose three simple problems in geochemical modeling that have nonunique solutions.

Changing the Basis

To this point we have assumed the existence of a basis of chemical components that corresponds to the system to be modeled. The basis, as discussed in the previous chapter, includes water, each mineral in the equilibrium system, each gas at known fugacity, and certain aqueous species. The basis serves two purposes: each chemical reaction considered in the model is written in terms of the members of the basis set, and the system’s bulk composition is expressed in terms of the components in the basis. Since we could not possibly store each possible variation on the basis, it is important for us to be able at any point in the calculation to adapt the basis to match the current system. It may be necessary to change the basis (make a basis swap, in modeling vernacular) for several reasons. This chapter describes how basis swaps can be accomplished in a computer model, and Chapter 9 shows how this technique can be applied to automatically balance chemical reactions and calculate equilibrium constants. The modeler first encounters basis swapping in setting up a model, when it may be necessary to swap the basis to constrain the calculation. The thermodynamic dataset contains reactions written in terms of a preset basis that includes water and certain aqueous species (Na+, Ca++, K+, Cl-, HCO-3, SO4- -, H+, and so on) normally encountered in a chemical analysis. Some of the members of the original basis are likely to be appropriate for a calculation. When a mineral appears at equilibrium or a gas at known fugacity appears as a constraint, however, the modeler needs to swap the mineral or gas in question into the basis in place of one of these species. Over the course of a reaction model, a mineral may dissolve away completely or become supersaturated and precipitate. In either case, the modeling software must alter the basis to match the new mineral assemblage before continuing the calculation. Finally, the basis sometimes must be changed in response to numerical considerations (e.g., Coudrain-Ribstein and Jamet, 1989). Depending on the numerical technique employed, the model may have trouble converging to a solution for the governing equations when one of the basis species occurs at small concentration.

Petroleum Reservoirs

In efforts to increase and extend production from oil and gas fields, as well as to keep wells operational, petroleum engineers pump a wide variety of fluids into the subsurface. Fluids are injected into petroleum reservoirs for a number of purposes, including: • Waterflooding, where an available fresh or saline water is injected into the reservoir to displace oil toward producing wells. • Improved Oil Recovery (IOR), where a range of more exotic fluids such as steam (hot water), caustic solutions, carbon dioxide, foams, polymers, surfactants, and so on are injected to improve recovery beyond what might be obtained by waterflooding alone. • Near-well treatments, in which chemicals are injected into producing and sometimes injector wells, where they are intended to react with the reservoir rock. Well stimulation techniques such as acidization, for example, are intended to increase the formation's permeability. Alternatively, producing wells may receive “squeeze treatments” in which a mineral scale inhibitor is injected into the formation. In this case, the treatment is designed so that the inhibitor sorbs onto mineral surfaces, where it can gradually desorb into the formation water during production. • Pressure management, where fluid is injected into oil fields in order to maintain adequate fluid pressure in reservoir rocks. Calcium carbonate may precipitate as mineral scale, for example, if pressure is allowed to deteriorate, especially in fields where formation fluids are rich in Ca++ and HCO3- and CO2 fugacity is high. In each of these procedures, the injected fluid can be expected to be far from equilibrium with sediments and formation waters. As such, it is likely to react extensively once it enters the formation, causing some minerals to dissolve and others to precipitate. Hutcheon (1984) appropriately refers to this process as “artificial diagenesis,” drawing an analogy to the role of groundwater flow in the diagenesis of natural sediments (see Chapter 19). Further reaction is likely if the injected fluid breaks through to producing wells and mixes there with formation waters. There is considerable potential, therefore, for mineral scale, such as barium sulfate (see the next section), to form during these procedures.

Stable Isotopes

Stable isotopes serve as naturally occurring tracers that can provide much information about how chemical reactions proceed in nature, such as which reactants are consumed and at what temperatures reactions occur. The stable isotopes of several of the lighter elements are sufficiently abundant and fractionate strongly enough to be of special usefulness. Foremost in importance are hydrogen, carbon, oxygen, and sulfur. The strong conceptual link between stable isotopes and chemical reaction makes it possible to integrate isotope fractionation into reaction modeling, allowing us to predict not only the mineralogical and chemical consequences of a reaction process, but also the isotopic compositions of the reaction products. By tracing the distribution of isotopes in our calculations, we can better test our reaction models against observation and perhaps better understand how isotopes fractionate in nature. Bowers and Taylor (1985) were the first to incorporate isotope fractionation into a reaction model. They used a modified version of EQ3/EQ6 (Wolery, 1979) to study the convection of hydrothermal fluids through the oceanic crust, along midocean ridges. Their calculation method is based on evaluating mass balance equations, as described in this chapter. As originally derived, however, the mass balance model has an important (and well acknowledged) limitation: implicit in its formulation is the assumption that fluid and minerals in the modeled system remain in isotopic equilibrium over the reaction path. This assumption is equivalent to assuming that isotope exchange between fluid and minerals occurs rapidly enough to maintain equilibrium compositions. We know, however, that isotope exchange in nature tends to be a slow process, especially at low temperature (e.g., O’Neil, 1987). This knowledge comes from experimental study (e.g., Cole and Ohmoto, 1986) as well as from the simple observation that, unless they have reacted together, groundwaters and minerals are seldom observed to be in isotopic equilibrium with each other. In fact, if exchange were a rapid process, it would be very difficult to interpret the origin of geologic materials from their isotopic compositions: the information would literally diffuse away. Lee and Bethke (1996) presented an alternative technique, also based on mass balance equations, in which the reaction modeler can segregate minerals from isotopic exchange. By segregating the minerals, the model traces the effects of the isotope fractionation that would result from dissolution and precipitation reactions alone.

Geochemical Buffers

Buffers are reactions that at least temporarily resist change to some aspect of a Fluid’s chemistry. A pH buffer, for example, holds pH to an approximately constant value, opposing processes that would otherwise drive the solution acid or alkaline.

Polythermal, Fixed, and Sliding Paths

In this chapter we consider how to construct reaction models that are somewhat more sophisticated than those discussed in the previous chapter, including reaction paths over which temperature varies and those in which species activities and gas fugacities are buffered. The latter cases involve the transfer of mass between the equilibrium system and an external buffer. Mass transfer in these cases occurs at rates implicit in solving the governing equations, rather than at rates set explicitly by the modeler. In Chapter 14 we consider the use of kinetic rate laws, a final method for defining mass transfer in reaction models. Polythermal reactions paths are those in which temperature varies as a function of reaction progress, ξ. In the simplest case, the modeler prescribes the temperatures T0 and Tf at the beginning and end of the reaction path. The model then varies temperature linearly with reaction progress. This type of model is sometimes called a “sliding temperature” path. The calculation procedure for a sliding temperature path is straightforward. In taking a reaction step, the model evaluates the temperature to be attained at the step’s end.

Mass Transfer

In previous chapters we have discussed the nature of the equilibrium state in geochemical systems: how we can define it mathematically, what numerical methods we can use to solve for it, and what it means conceptually. With this chapter we begin to consider questions of process rather than state. How does a fluid respond to changes in composition as minerals dissolve into it, or as it mixes with other fluids? How does a fluid evolve in response to changing temperature or variations in the fugacity of a coexisting gas? In short, we begin to consider reaction modeling. In this chapter we consider how to construct reactions paths that account for the effects of simple reactants, a name given to reactants that are added to or removed from a system at constant rates. We take on other types of mass transfer in later chapters. Chapter 12 treats the mass transfer implicit in setting a species’ activity or gas’ fugacity over a reaction path. In Chapter 14 we develop reaction models in which the rates of mineral precipitation and dissolution are governed by kinetic rate laws. Simple reactants are those added to (or removed from) the system at constant rates over the reaction path. As noted in Chapter 2, we commonly refer to such a path as a titration model, because at each step in the process, much like in a laboratory titration, the model adds an aliquot of reactant mass to the system. Each reactant Ar is added at a rate nr, expressed in moles per unit reaction progress, ξ. Negative values of nr, of course, describe the removal rather than the addition of the reactant. Since ξ is unitless and varies from zero at the start of the path to one at the end, we can just as well think of nr as the number of moles of the reactant to be added over the reaction path. A simple reactant may be an aqueous species (including water), a mineral, a gas, or any entity of known composition.

Surface Complexation

An important consideration in constructing certain types of geochemical models, especially those applied to environmental problems, is to account for the sorption of ions from solution onto mineral surfaces. Metal oxides and aluminosilicate minerals, as well as other phases, can sorb electrolytes strongly because of their high reactivities and large surface areas (e.g., Davis and Kent, 1990). When a fluid comes in contact with minerals such as iron or aluminum oxides and zeolites, sorption may significantly diminish the mobility of dissolved components in solution, especially those present in minor amounts. Sorption, for example, may retard the spread of radionuclides near a radioactive waste repository or the migration of contaminants away from a polluting landfill. In acid mine drainages, ferric oxide sorbs heavy metals from surface water, helping limit their downstream movement (see Chapter 23). A geochemical model useful in investigating such cases must provide an accurate assessment of the effects of surface reactions. Many of the sorption theories now in use are too simplistic to be incorporated into a geochemical model intended for general use. To be useful in modeling electrolyte sorption, a theory must account for the electrical charge on the mineral surface and provide for mass balance on the sorbing sites. In addition, an internally consistent and sufficiently broad database of sorption reactions must accompany the theory. The Freundlich and Langmuir theories, which use distribution coefficients Kd to set the ratios of sorbed to dissolved ions, are applied widely in groundwater studies (Domenico and Schwartz, 1990) and used with considerable success to describe sorption of uncharged organic molecules (Adamson, 1976). The models, however, do not account for the electrical state of the surface, which varies sharply with pH, ionic strength, and solution composition. Freundlich theory prescribes no concept of mass balance, so that a surface might be predicted to sorb from solution without limit. Both theories require that distribution coefficients be determined experimentally for individual fluid and rock compositions, and hence both theories lack generality. Ion exchange theory (Stumm and Morgan, 1981; Sposito, 1989) suffers from similar limitations. Surface complexation models, on the other hand, account explicitly for the electrical state of the sorbing surface (e.g., Adamson, 1976; Stumm, 1992).

Solving for the Equilibrium State

In Chapter 3, we developed equations that govern the equilibrium state of an aqueous fluid and coexisting minerals. The principal unknowns in these equations are the mass of water nw, the concentrations mi of the basis species, and the mole numbers nk of the minerals. If the governing equations were linear in these unknowns, we could solve them directly using linear algebra. However, some of the unknowns in these equations appear raised to exponents and multiplied by each other, so the equations are nonlinear. Chemists have devised a number of numerical methods to solve such equations (e.g., van Zeggeren and Storey, 1970; Smith and Missen, 1982). All the techniques are iterative and, except for the simplest chemical systems, require a computer. The methods include optimization by steepest descent (White et al., 1958; Boynton, 1960) and gradient descent (White, 1967), back substitution (Kharaka and Barnes, 1973; Truesdell and Jones, 1974), and progressive narrowing of the range of the values allowed for each variable (the monotone sequence method; Wolery and Walters, 1975). Geochemists, however, seem to have reached a consensus (e.g., Karpov and Kaz’min, 1972; Morel and Morgan, 1972; Crerar, 1975; Reed, 1982; Wolery, 1983) that Newton-Raphson iteration is the most powerful and reliable approach, especially in systems where mass is distributed over minerals as well as dissolved species. In this chapter, we consider the special difficulties posed by the nonlinear forms of the governing equations and discuss how the Newton-Raphson method can be used in geochemical modeling to solve the equations rapidly and reliably. The governing equations are composed of two parts: mass balance equations that require mass to be conserved, and mass action equations that prescribe chemical equilibrium among species and minerals.

Acid Drainage

Acid drainage is a persistent environmental problem in many mineralized areas. The problem is especially pronounced in areas that host or have hosted mining activity (e.g., Lind and Hem, 1993), but it also occurs naturally in unmined areas. The acid drainage results from the weathering of sulfide minerals that oxidize to produce hydrogen ions and contribute dissolved metals to solution. These acidic waters are toxic to plant and animal life, including fish and aquatic insects. Streams affected by acid drainage may be rendered nearly lifeless, their stream beds coated with unsightly yellow and red precipitates of oxy-hydroxide minerals. In some cases, the heavy metals in acid drainage threaten water supplies and irrigation projects. Where acid drainage is well developed and extensive, the costs of remediation can be high. In the Summitville, Colorado district (USA), for example, efforts to limit the contamination of fertile irrigated farmlands in the nearby San Luis Valley and protect aquatic life in the Alamosa River will cost an estimated $100 million or more (Plumlee, 1994a). Not all mine drainage, however, is acidic or rich in dissolved metals (e.g., Ficklin et al., 1992; Mayo et al., 1992; Plumlee et al., 1992). Drainage from mining districts in the Colorado Mineral Belt ranges in pH from 1.7 to greater than 8 and contains total metal concentrations ranging from as low as about 0.1 mg/kg to more than 1000 mg/kg. The primary controls on drainage pH and metal content seem to be (1) the exposure of sulfide minerals to weathering, (2) the availability of atmospheric oxygen, and (3) the ability of nonsulfide minerals to buffer acidity. In this chapter we construct geochemical models to consider how the availability of oxygen and the buffering of host rocks affect the pH and composition of acid drainage. We then look at processes that can attenuate the dissolved metal content of drainage waters. Acid drainage results from the reaction of sulfide minerals with oxygen in the presence of water. As we show in this section, water in the absence of a supply of oxygen gas becomes saturated with respect to a sulfide mineral after only a small amount of the mineral has dissolved.