Soil Science Society of America Journal 66:474-478 (2002)
© 2002 Soil Science Society of America
DIVISION S-2—NOTES
Evaluation of numerical techniques applied to soil solution speciation including cation exchange
Peter J. Vaughan*
Salinity Laboratory USDA-ARS, 450 W. Big Springs Rd., Riverside, CA 92507
* Corresponding author (pvaughan{at}ussl.ars.usda.gov)
 |
ABSTRACT
|
|---|
Most existing models of soil solution speciation utilize either Newton–Raphson or Picard iteration to obtain a numerical solution to the nonlinear set of algebraic equations expressing mole balance and charge balance of free ions for major species as well as cation exchange. A computer program was written to test speed and accuracy of these and other methods. Picard iteration was fastest but produced a mean relative error (MRE) for mole and charge balance of 0.03 with no further convergence after a few iterations. A tensor (quadratic) method and two simplex methods converged to the correct result. The tensor method was preferable because the final result was more accurate; the rate of convergence was faster by 10 to 100 times, and convergence occurred for all compositions tested up to an ionic strength of 0.25. These results point out the value of testing various algorithms prior to implementation of speciation code applied to soils.
Abbreviations: CEC, cation-exchange capacity • MRE, mean relative error • RMS, root mean squares
|
INTRODUCTION
|
|---|
NUMERICAL CALCULATION of equilibrium chemical speciation is a standard procedure that has been the goal of many different models including WATEQ4F (Ball and Nordstrom, 1991), PhreeqeC (Parkhurst and Appelo, 1999), and EQ3NR (Wolery, 1992). These models include other types of chemical reactions in addition to speciation within the aqueous phase such as cation exchange, dissolution and precipitation of solid phases, and reactions between surface species and the bulk solution. For the bulk solution, a set of master species is normally constructed. These comprise a minimum set of species from which all other species in the system can be obtained by reaction. A numerical solution to the mole balance and charge balance equations is commonly obtained by successive approximation of the concentration vector of master species by the Newton–Raphson method. This method relies on calculation of the Jacobian, a matrix of the partial derivatives of the vector of residual functions with respect to the master species concentrations. The Newton–Raphson method is known to have poor convergence or even failure to converge for certain initial guesses of the concentration vector (Parkhurst and Appelo, 1999; Schnabel and Frank, 1984). Furthermore, the convergence criteria for the sequence of successive approximations, by Newton–Raphson, to the correct solution are not necessarily met unless the vector for the initial guess is within a specified radius of the correct solution (Holstad, 1999). Given the difficulties associated with the Newton–Raphson method as applied to speciation calculations, it appears worthwhile to examine the effectiveness of various numerical algorithms applied to the solution of speciation problems common in soil science. Effectiveness needs to be assessed both in terms of the accuracy of the numerical solution and its speed. The importance of accuracy is obvious but speed requires some clarification. A common soil science application of speciation codes would be solute transport problems as treated by a finite-element model. In these models speed is critical because of the large fraction of total execution time that is spent computing the equilibrium speciation. For example, profiling tests of the Unsatchem multicomponent transport model showed this fraction commonly exceeding 0.8 (P. Vaughan, unpublished data, 2001).
Solution speciation problems are normally formulated as a set of nonlinear equations representing relationships among master and secondary species (Parkhurst, 1997; Morel and Morgan, 1972). The relationships among all species are expressed through consideration of mole balance, electroneutrality, and the mass action laws representing reactions. For soil solutions, cation exchange reactions also need to be considered. This set of mostly nonlinear equations can be solved numerically by minimizing a residual function that provides quantification of the total error associated with each successive approximation.
The Newton–Raphson and Levenberg–Marquardt methods rely on the Jacobian to make an improved estimate of the unknown (Press et al., 1986). Some other numerical techniques can also address the speciation problem but do not require direct computation of the Jacobian. These include Picard iteration, the Nelder–Mead simplex algorithm (Nelder and Mead, 1965), global methods such as simulated annealing techniques, and a tensor method that utilizes a quadratic expression to compute successive iterates (Schnabel and Frank, 1984). Several methods based on these various algorithms were tested on a problem of solving for master species' concentrations when both cation exchange and solution speciation are considered.
|
Materials and Methods
|
|---|
The problem to be solved consists of a set of mole balance equations for each of the master species. These include master species appearing as free ions in solution, components in secondary species and on exchange sites on the surfaces of solids. An additional constraint is that of electroneutrality in the bulk solution.
Cation-exchange reactions can be written in various ways including the Gapon, Vanselow, and Gaines-Thomas formulations (Sposito, 1981). This paper utilizes the Gapon formulation in which the reaction is represented in terms of equivalents.
 | [1] |
where X represents the concentration of either the m or n cation on the exchange phase (mmolc) and M, N represent the activities of each cation in solution (Robbins et al., 1980). The mass action expressions, of the form
 | [2] |
can be combined with an equation expressing the cation-exchange capacity (CEC) as the sum of the exchangeable cations (mmolc) to obtain expressions for the exchange concentrations (Robbins et al., 1980),
In this equation the three Gapon selectivity coefficients were defined so that they appeared in either numerator or denominator of the various terms. From the standpoint of convenience it's easier to define the coefficients so that one exchange species appears consistently in the numerator of the mass action law. Choosing X1/2Ca gives the following expression for Ca-Mg exchange
 | [4] |
that is the reciprocal of the expression for K1 given by Robbins et al. (1980). Expressions for the Gapon selectivity coefficients for Ca-Na and Ca-K exchange are identical to those of Robbins et al. (1980). The exchange-phase concentration can be written as
 | [5] |
where (wj) is the activity of the jth cation, aj is the stoichiometric coefficient and Kj is the selectivity coefficient. For j= 1 the selectivity coefficient represents the exchange of 1/2Ca2+ with X1/2Ca and is equal to one by definition.
Mass action expressions for ion pairs in solution can be represented by:
 | [6] |
where (wi), (yi), (ci) are cation, anion, and ion pair activities, respectively. The superscripts j and k are stoichiometric coefficients and the Ki is an equilibrium constant. Carbonate chemistry is computed on the assumption that pCO2 is externally fixed and that total alkalinity
 | [7] |
is conserved (Simunek et al., 1996). The brackets in Eq [7] signify molality.
The residuals for the mole balance equations for each of the master species and overall electroneutrality form a residual error vector that must be minimized to obtain the vector of master species concentrations. For the algorithms discussed here a single-valued objective function was required; therefore, the root mean square (RMS) of the residual error vector was computed as this value. Because the exchangeable concentrations are expressed as mmolc/kg soil, conversion of the exchangeable concentration to a hypothetical molality is required before summation with the molality of the corresponding species in solution,
 | [8] |
In this expression, 
is the soil bulk density (kg m-3), 
is the volumetric water content, zi is the ionic charge, and Xi is the exchangeable-cation concentration (mmolc kg-1 soil).
Activity Coefficients
Activity coefficients for species in solution were computed from ionic strength using the extended Debye-Huckel approximation (Truesdell and Jones, 1974).
 | [9] |
The parameters a and b are species-specific whereas A and B are dependent only on the dielectric constant, temperature, and solution density. The variable, zi, denotes the ionic charge of each species and I is the ionic strength of the solution (mol kg-1).
Data Sets
Five synthetic data sets were created to provide a range of saturation paste-extract compositions and CECs typical of soils (Table 1). The ionic strength of the solution should be <0.2 to justify utilization of the extended Debye-Huckel approximation for computation of the activity coefficients using specific ion coefficients (Truesdell and Jones, 1974). The last solution listed in Table 1 had an ionic strength of 0.247 that was greater than the recommended range. However, the objective of this exercise was testing a range of possible compositions in evaluating the performance and stability of the various algorithms. The Gapon selectivity coefficients, as defined in Eq. [5], were 
.
Numerical Tests
Performance of the various methods was based on a combination of speed and accuracy. The speed component was measured by determination of the central processing unit (CPU) time required from start to completion of the iterative portion of the program. Assessment of accuracy was based on relative error.
All algorithms were coded in Fortran 77 and all runs were performed on a 50 Mhz Sun SPARC 20 1(Sun MicroSystems, Palo Alto, CA). Five methods tested included two implementations of the Nelder–Mead simplex algorithm (a modified simplex algorithm - Cobyla2), a standard Nelder–Mead implementation, a global solver (Toms667), a tensor method (Tensolve) and a Picard (fixed point) iteration. The subroutine for each method was passed a starting composition including total free concentrations of master species, zero concentration for secondary species and cation concentrations computed for a single exchanger. The initial exchangeable concentrations were calculated from Eq. [5] assuming that the free concentrations of cations were the totals for the master cations. The final RMS of the mole and charge balance error was computed by a separate subroutine that was identical for all five methods.
The elapsed time allotted for a single test was controlled to examine the tradeoff between elapsed time and accuracy. This control was implemented in various ways for the different methods. For example, control of the simplex methods was accomplished by adjusting the maximum number of function evaluations while the tensor method was controlled through adjustment of the tolerance for the RMS of the combined mole and charge balance errors.
To ensure that numerical results were correct, the concentrations of the master and secondary species were used to back calculate the equilibrium constants. Also, the Gapon selectivity coefficients and CEC were back calculated from the master cations and exchangeable-cation concentrations.
|
Results and Discussion
|
|---|
Several algorithms were studied using the MatLab program on a desktop personal computer. The Levenberg–Marquardt and Gauss–Newton methods were programmed to include a direct calculation of the Jacobian. A comparison of the performance of these algorithms when the Jacobian was directly calculated and when it was approximated by finite differences indicated, surprisingly, that the approximation of the Jacobian resulted in faster convergence with less likelihood of the solution becoming trapped in a false minimum. Direct computation of the Jacobian also had increased the likelihood of an unsuccessful result given the same choice of initial conditions. For these reasons the algorithms requiring computation of the Jacobian were dropped from further consideration.
The remaining five algorithms were tested on the workstation. Results of the first test demonstrated that the toms667 global solver was not an appropriate choice for this problem because of its slow convergence (Fig. 1)
. Subsequent tests were performed only on the other four methods.
Among the remaining methods there were large differences in both accuracy and speed of the numerical solution. The timings reported here were for Composition 2 with CEC = 100 mmolc kg-1 soil (Table 1). The fastest algorithm was the Picard iteration with an elapsed time of 1.4 ms for the first iteration. When solution speciation was calculated with no cation exchange the Picard iteration provided a fast and accurate solution. Inclusion of the cation exchange, however, resulted in lack of convergence with a minimum MRE in mole balance and charge balance of 
0.03 for Composition 2 (Table 1). The continuation of Picard iteration did not provide further improvement in the accuracy. For the remaining compositions, as ionic strength increased above 0.123 and the CEC was also increased, the Picard method showed no convergence. It should be noted that Picard iteration was done in stages with the solution speciation alternating with cation exchange and pH calculations. It is certainly possible that some rearrangement of this code could potentially produce better results. This does not apply to the other methods, which are all simultaneous solutions.
The two simplex methods converged slowly but numerical solutions with MRE of 0.001 or better were obtained for Composition 2 (Fig. 1). A typical comparison of elapsed time to achieve the same MRE of 1.7 x 10-5 was (tensor method, 0.173 s; Cobyla2, 7.6 s; Nelder-Mead, 40 s). The tensor method showed rapid convergence in all tests and consistently achieved a MRE of better than 1 x 10-12. Results obtained for the other four solution compositions were similar to those obtained for Composition 2. A significant difference in MRE among the different solution compositions was a progressively lower accuracy obtained by the standard Nelder–Mead method with increasing ionic strength. The modified simplex method (Cobyla2) did not suffer from this problem. The correct back calculation of equilibrium constants, Gapon selectivity coefficients, and CEC from master species concentrations obtained from the tensor method demonstrated that the correct numerical solution of the equations had been obtained.
The relative performance of the methods can be evaluated by comparing the MRE at some fixed elapsed time for Composition 2 (Table 2). The two simplex methods actually have higher MRE than the starting guess and substantially greater error than the Picard iteration at 1.0 s. At 1 s elapsed time, the tensor method had obtained a MRE of 2 x 10-12. At 10 s elapsed time, the Cobyla2 simplex method obtained a MRE of 1.1 x 10-4, a significant improvement over the Picard iteration. The standard Nelder–Mead algorithm, however, was less accurate than Picard at both 1 and 10 s.
These results indicate that the choice of algorithm depends on the desired accuracy and speed. If only a rough approximation is needed, ionic strength is <0.2, and speed is a critical factor then the Picard iteration would be the logical choice. The simplex methods should only be used when their convergence can be assured by allocating sufficient computation time. The Cobyla2 method performed consistently better than standard Nelder–Mead. Both of these methods have a weakness of slow convergence early in the process. Overall, the tensor method performed best among the algorithms tested.
|
Conclusions
|
|---|
Various numerical techniques were tested for suitability in the problem of equilibrium speciation when cation exchange is also included. Given this problem and typical starting total concentrations for soil paste extracts, the methods requiring computation of the Jacobian such as Newton–Raphson and Levenberg–Marquardt were found to have convergence that was highly dependent on the starting guess. Methods that consistently obtained accurate solutions included two simplex methods and the tensor method. Picard iteration provided the fastest computation of a rough approximation having a relative error in mean mole balance of 0.03. The two simplex methods converged more slowly than the tensor method by a factor of 10 to 1000 when obtaining the same accuracy. The tensor method was judged to be the best choice considering both speed and accuracy.
|
NOTES
|
|---|
Sponsoring Organization: Agricultural Research Service, USDA.
1 The use of brand names in this report is for identification purposes only and does not constitute endorsement by the USDA. 
Received for publication May 21, 2001.
|
REFERENCES
|
|---|
- Ball, J.W., and D.K Nordstrom. 1991. User's manual for WATEQ4F, with revised thermodynamic data base and test cases for calculating speciation of major, trace, and redox elements in natural waters. U.S. Geological Survey Report 91-183. [online] Available at: http://water.usgs.gov/software/wateq4f.html (verified 19 Nov. 2001).
- Holstad, A. 1999. Numerical solution of nonlinear equations in chemical speciation calculations. Comput. Geosci. 3:229–257.
- Morel, F., and J. Morgan. 1972. A numerical method for computing equilibria in aqueous chemical systems. Environ. Sci. Technol. 6:58–67.
- Nelder, J.A., and R. Mead. 1965. A simplex method for function minimization. Comput. J. 7:308–313.
- Parkhurst, D.L. 1997. Geochemical mole balance modeling with uncertain data. Water Resour. Res. 33:1957–1970.
- Parkhurst, D.L. and C.A.J. Appelo. 1999. User's guide to PhreeqC (Version 2)—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Water-Resources Investigations Rep. 99-4259, U.S. Geological Survey, Denver, CO.
- Press, W.H., Flannery, B.P., Teukolsky, S.A., and W.T. Vetterling. 1986. Numerical Recipes. Cambridge University Press, Cambridge.
- Robbins, C.W., Jurinak, J.J., and R.J. Wagenet. 1980. Calculating cation exchange in a salt transport model. Soil Sci. Soc. Am. J. 44:1195–1200.[Abstract/Free Full Text]
- Schnabel, R.B., and P.D. Frank. 1984. Tensor methods for nonlinear equations. SIAM J. Num. Anal. 21:815–843.
- Simunek, J., Suarez, D.L., and M. Sejna. 1996. The UNSATCHEM software package for simulating the one-dimensional variably saturated water flow, heat transport, carbon dioxide production and transport, and multicomponent solute transport with major ion equilibrium and kinetic chemistry, Ver. 2.0. Research Rep. 141. George E. Brown, Jr. Salinity Laboratory USDA-ARS, Riverside, CA.
- Sposito, G. 1981. The thermodynamics of soil solutions. Clarendon, Oxford.
- Truesdell, A.H., and B.F. Jones. 1974. WATEQ, A computer program for calculating chemical equilibria in natural waters. J. Res. U.S. Geol. Surv. 2:233–248.
- Wolery, T.J. 1992. EQ3NR, A computer program for geochemical aqueous speciation-solubility calculations: Theoretical manual, user's guide, and documentation, Version 7.0. Lawrence-Livermore National Laboratory Report UCRL-MA-110662 PT III, Livermore, CA.
TOP
ABSTRACT
INTRODUCTION
) vs. log10(t) where