A typical environmental study usually have several components and need to solve dozens of conservation equations to convergence tolerance of a very small value due to the low concentration of pollutants, which will take the computer tremendous time. A new philosophy of solving the conservation equations is employed in this study. The new method only solves the conservation of minor components to the strict convergence tolerance due to their small concentrations. When we solve the equations of main components, we solve the equations to a relaxed convergence tolerance since their concentrations are much larger than those of minor components.

The new method was implemented and tested in a commercial environmental simulator that includes speciesdependent convergence tolerance, such that all conservations are solved to similar accuracy. Various field test cases are discussed; whose run times are faster than the original one with excellent comparison in pollutants distribution. The new method of solving the conservation equations greatly reduced the cost of environmental studies and made environmental simulation studies more manageable.

Any fluid flow model in multi-component geothermal simulators is based on the conservation of mass, energy and momentum. Here we take the widely used, commercial multicomponent geothermal simulator TETRAD as an example. The conservation of component i (for i=1, total number of components, or Nc) can be written as:

Equation (1) (Available in full paper)

Equation (2) (Available in full paper)

The conservation equations of energy and momentum are very similar to the above equation (1). The first term in equation (1) above is the accumulation of component i, the second term is the net flux (which includes advection and dispersion), qi is the source/sink of component i, and Ri is formation/destruction of i from reaction or decay. Individual terms are defined below: ωij are the mole fraction of component i in phase j (for each of j=1, Np phases, s=solid) ρj are phase densities Sj are phase saturations φis porosity uj are Darcy velocities of phase j K_{ij} is the dispersion tensor.

A general discussion of these equations and relevant constitutive equations is given by Vinsome and Shook (1993) and Shook (1995).

The conservation equation given above (along with all constitutive relations) is solved for each chemical component at each grid node in the domain of interest. As an example, Shook and Faulder (1991) showed that a total of 27 equations are required at each grid node for a two-component geothermal (i.e., not isothermal) problem. Equations at each grid node are coupled through standard finite differencing; thus the equation set is an N-by-N, banded, block matrix with Nc-by-Nc blocks for each of N grid-blocks.

The solution of these Nc coupled equations is accomplished through a Newton-Raphson procedure as shown below.

Equation (3) (Available in full paper)

where J (x) is the Jacobian matrix (partial derivatives of the conservation equations with respect to primary variables) that is implicit in time. That is, the coefficients in the conservation equations Ei are functions of primary variables at t+ ⊗t, are therefore currently unknown.