Practical aspects of three-dimensional modeling of groundwater flow in heterogeneous aquifer systems are investigated using a finite element approach. Particular attention is given to the properties of the conductance matrix and the efficiency of the conjugate gradient method with different preconditioners: diagonal scaling, incomplete Cholesky decomposition, incomplete factorization, and modified incomplete factorization. It is shown that for hexahedral trilinear finite elements the resulting matrix is, except for cube-shaped elements, never diagonally dominant, which restricts the existence of several preconditioners. Numerical comparison of several test problems, including hypothetical and field applications with different degrees of heterogeneity, show that the incomplete Cholesky and the incomplete factorization preconditioners, if they exist, are more efficient than diagonal scaling with respect to both rate of convergence and overall computing time, but diagonal scaling can be considered superior because it is always possible. An M matrix transformation is proposed which guarantees the existence of all preconditioners. Numerical comparison of the test problems shows that this technique is very effective. From the resulting preconditioners, the incomplete Cholesky and the incomplete factorization are shown to be the most efficient, but the latter is superior from the point of view of computer storage and is recommended for all practical applications. ¿American Geophysical Union 1994