Geochemists have long recognized the need for computational models to trace the progress of reaction processes, both natural and artificial. Given a process involving many individual reactions (possibly thousands), some of which yield products that provide reactants for others, how can we know which reactions are important, how far each will progress, what overall reaction path will be followed, and what the path's endpoint will be?

These questions can be answered reliably by hand calculation only in simple cases. Geochemists are increasingly likely to turn to quantitative modeling techniques to make their evaluations, confirm their intuitions, and spark their imaginations.

Computers were first used to solve geochemical models in the 1960s, but the new modeling techniques disseminated rather slowly through the practice of geochemistry. Even today, many geochemists consider modeling to be a “black art,” perhaps practiced by digital priests muttering mantras like “Newton–Raphson” and “Runge–Kutta” as they sit before their cathode ray altars. Others show little fear in constructing models but present results in a way that adds little understanding of the problem considered. Someone once told me, “Well, that's what came out of the computer!”

A large body of existing literature describes either the formalism of numerical methods in geochemical modeling or individual modeling applications. Few references, however, provide a perspective of the modeling specialty, and some that do are so terse and technical as to discourage the average geochemist. Hence, there are few resources to which someone wishing to construct a model without investing a career can turn.

Conveniently, perhaps even miraculously, the equations developed in Chapter 5 to accomplish basis swaps can be used to balance chemical reactions automatically. Once the equations have been coded into a computer program, there is no need to balance reactions, compute equilibrium constants, or even determine equilibrium equations by hand. Instead, these procedures can be performed quickly and reliably on a small computer.

To balance a reaction, we first choose a species to appear on the reaction's left side, and express that species' composition in terms of a basis B. The basis might be a list of the elements in the species' stoichiometry, or an arbitrary list of species that combine to form the left-side species. Then we form a second basis B′ composed of species that we want to appear on the reaction's right side. To balance the reaction, we calculate the transformation matrix relating basis B′ to B, following the procedures in Chapter 5. The transformation matrix, in turn, gives the balanced reaction and its equilibrium constant.

Calculation procedure

Two methods of balancing reactions are of interest. We can balance reactions in terms of the stoichiometries of the species considered. In this case, the existing basis B is a list of elements and, if charged species are involved, the electron e−. Alternatively, we may use a dataset of balanced reactions, such as the llnl database. Basis B, in this case, is the one used in the database to write reactions. We will consider each possibility in turn.

Many times in geochemical modeling we want to understand not only what reactions proceed in an open chemical system, but where they take place (e.g., Steefel et al., 2005). In a problem of groundwater contamination, for example, we may wish to know not only the extent to which a contaminant might sorb, precipitate, or degrade, but how far it will migrate before doing so.

To this point, we have discussed how to model the reactions occurring within a single representative volume, as shown in Figure 2.1. Such configurations are sometimes called lumped parameter models, because the properties of the entire system are represented by a single set of values. There is one pH for the entire volume, one ionic strength, a single mass for each chemical component, and so on. In studying reaction within flowing groundwater, in contrast, we may want to build a distributed parameter model, a model in which the properties vary across the system.

To construct models of this sort, we combine reaction analysis with transport modeling, the description of the movement of chemical species within flowing groundwater, as discussed in the previous chapter (Chapter 20). The combination is known as reactive transport modeling, or, in contaminant hydrology, fate and transport modeling.

Mathematical model

A reactive transport model, as the name implies, is reaction modeling implemented within a transport simulation. It may be thought of as a reaction model distributed over a groundwater flow. In other words, we seek to trace the chemical reactions that occur at each point in space, accounting for the movement of reactants to that point, and reaction products away from it.

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 temperature 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).

Hydrothermal fluids, hot groundwaters that circulate within the Earth's crust, play central roles in many geological processes, including the genesis of a broad variety of ore deposits, the chemical alteration of rocks and sediments, and the origin of hot springs and geothermal fields. Many studies have been devoted to modeling how hydrothermal fluids react chemically as they encounter wall rocks, cool, boil, and mix with other fluids. Such modeling proliferated in part because hydrothermal fluids are highly reactive and because the reaction products are commonly well preserved, readily studied, and likely to be of economic interest. Further impetus was provided by the development of reliable modeling software in the 1970s, a period of concern over the availability of strategic and critical minerals and of heightened interest in economic geology and the exploitation of geothermal energy.

As a result, many of the earliest and most imaginative applications of geochemical modeling, beginning with Helgeson's (1970) simulation of ore deposition in hydrothermal veins and the alteration of nearby country rock, have addressed the reaction of hydrothermal fluids. For example, Reed (1977) considered the origin of a precious metal district; Garven and Freeze (1984), Sverjensky (1984, 1987), and Anderson and Garven (1987) studied the role of sedimentary brines in forming Mississippi Valley-type and other ore deposits; Wolery (1978), Janecky and Seyfried (1984), Bowers et al. (1985), and Janecky and Shanks (1988) simulated hydrothermal interactions along the midocean ridges; and Drummond and Ohmoto (1985) and Spycher and Reed (1988) modeled how fluid boiling is related to ore deposition.

Diagenesis is the set of processes by which sediments evolve after they are deposited and begin to be buried. Diagenesis includes physical effects such as compaction and the deformation of grains in the sediment (or sedimentary rock), as well as chemical reactions such as the dissolution of grains and the precipitation of minerals to form cements in the sediment's pore space. The chemical aspects of diagenesis are of special interest here.

Formerly, geologists considered chemical diagenesis to be a process by which the minerals and pore fluid in a sediment reacted with each other in response to changes in temperature, pressure, and stress. As early as the 1960s and especially since the 1970s, however, geologists have recognized that many diagenetic reactions occur in systems open to groundwater flow and mass transfer. The reactions proceed in response to a supply of reactants introduced into the sediments by flowing groundwater, which also serves to remove reaction products.

Hay (1963, 1966), in studies of the origin of diagenetic zeolite, was perhaps the first to emphasize the effects of mass transport on sediment diagenesis. He showed that sediments open to groundwater flow followed reaction pathways different from those observed in sediments through which flow was restricted. Sibley and Blatt (1976) used cathodoluminescence microscopy to observe the Tuscarora orthoquartzite of the Appalachian basin. The almost nonporous Tuscarora had previously been taken as a classic example of pressure welding, but the microscopy demonstrated that the rock is not especially well compacted but, instead, tightly cemented. The rock consists of as much as 40% quartz (SiO2) cement that was apparently deposited by advecting groundwater.

To this point we have measured reaction progress parametrically in terms of the reaction progress variable ξ, which is dimensionless. When in Chapter 13 we reacted feldspar with water, for example, we tied reaction progress to the amount of feldspar that had reacted and expressed our results along that coordinate. Studying reactions in this way is in many cases perfectly acceptable. But what if we want to know how much time it took to reach a certain point along the reaction path? Or, when modeling the reaction of granite with rainwater, how can we set the relative rates at which the various minerals in the granite dissolve? In such cases, we need to incorporate reaction rate laws from the field of geochemical kinetics.

The differences between the study of thermodynamics and kinetics might be illustrated (e.g., Lasaga, 1981a) by the analogy of rainfall on a mountain. On the mountaintop, the rainwater contains a considerable amount of potential energy. With time, it flows downhill, losing energy (to be precise, losing hydraulic potential, the mechanical energy content of a unit mass of water; Hubbert, 1940), until it eventually reaches the ocean, its lowest possible energy level. The thermodynamic interpretation of the process is obvious: the water seeks to minimize its energy content.

But how long will it take for the rainfall to reach the ocean? The rain might enter a swift mountain stream, flow into a river, and soon reach the sea. It might infiltrate the subsurface and migrate slowly through deep aquifers until it discharges in a distant valley, thousands of years later.

Redox reactions in the geochemical environment, as discussed in previous chapters (Chapters 7 and 17), are commonly in disequilibrium at low temperature, their progress described by kinetic rate laws. The reactions may proceed in solution homogeneously or be catalyzed on the surface of minerals or organic matter. In a great many cases, however, they are promoted by the enzymes of the ambient microbial community.

In this chapter, we consider how the microbial community catalyzes redox reactions, perhaps changing in size and composition as it does. The kinetics of such reactions are of special interest, because of the close relationship between geochemical conditions and microbial ecology. The microbes promote reactions that change geochemical conditions, many times significantly, and the geochemistry controls the nature of the microbial community that can exist in a given environment.

From the geochemist's perspective, a microbe can be thought of as a selfreplicating bundle of enzymes. Microbes use their collections of enzymes to catalyze redox reactions, harvesting some of the energy released for their own purposes. Since microbial growth increases the ability of the community to catalyze redox reactions, and catalyzing the reactions provides microbes the energy they need to grow, a microbially mediated reaction is by nature autocatalytic.

Microbial respiration and fermentation

Chemosynthetic microorganisms derive the energy they need to live and grow from chemical species in their environment, reaping the benefits of the redox disequilibrium characteristic of geochemical environments (see Chapter 7).

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.

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 14 treats the mass transfer implicit in setting a species' activity or gas' fugacity over a reaction path. In Chapter 16 we develop reaction models in which the rates of mineral precipitation and dissolution are governed by kinetic rate laws.

Simple reactants

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, ξ.

As geochemists, we frequently need to describe the chemical states of natural waters, including how dissolved mass is distributed among aqueous species, and to understand how such waters will react with minerals, gases, and fluids of the Earth's crust and hydrosphere. We can readily undertake such tasks when they involve simple chemical systems, in which the relatively few reactions likely to occur can be anticipated through experience and evaluated by hand calculation. As we encounter more complex problems, we must rely increasingly on quantitative models of solution chemistry and irreversible reaction to find solutions.

The field of geochemical modeling has grown rapidly since the early 1960s, when the first attempt was made to predict by hand calculation the concentrations of dissolved species in seawater. Today's challenges might be addressed by using computer programs to trace many thousands of reactions in order, for example, to predict the solubility and mobility of forty or more elements in buried radioactive waste.

Geochemists now use quantitative models to understand sediment diagenesis and hydrothermal alteration, explore for ore deposits, determine which contaminants will migrate from mine tailings and toxic waste sites, predict scaling in geothermal wells and the outcome of steam-flooding oil reservoirs, solve kinetic rate equations, manage injection wells, evaluate laboratory experiments, and study acid rain, among many examples. Teachers let their students use these models to learn about geochemistry by experiment and experience.

Many hundreds of scholarly articles have been written on the modeling of geochemical systems, giving mathematical, geochemical, mineralogical, and practical perspectives on modeling techniques.