Accuracy of Galerkin finite elements for groundwater flow simulations in two and three-dimensional triangulations

In standard finite element simulations of groundwater flow the correspondence between hydraulic head gradients and groundwater fluxes is represented by the stiffness matrix. In two-dimensional problems the use of linear triangular elements on Delaunay triangulations guarantees a stiffness matrix of type M. This implies that the local numerical fluxes are physically consistent with Darcy's law. This condition is fundamental to avoid the occurrence of local maxima or minima, and is of crucial importance when the calculated flow field is used in contaminant transport simulations or pathline evaluation. In three spatial dimensions, the linear Galerkin approach on tetrahedra does not lead to M-matrices even on Delaunay meshes. By interpretation of the Galerkin approach as a subdomain collocation scheme, we develop a new approach (OSC, orthogonal subdomain collocation) that is shown to produce M-matrices in three-dimensional Delaunay triangulations. In case of heterogeneous and anisotropic coefficients, extra mesh properties required for M-stiffness matrices will also be discussed.