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.