Saltwater intrusion in aquifers: Development and testing of a three‐dimensional finite element model

Abstract

A three‐dimensional finite element model is developed for the simulation of saltwater intrusion in single and multiple coastal aquifer systems with either a confined or phreatic top aquifer. The model formulation is based on two governing equations, one for fluid flow and the other for salt transport. Density coupling of these equations is accounted for and handled using a Picard sequential solution algorithm with special provisions to enhance convergence of the iterative solution. Flexibility in the formulation allows for either three‐dimensional simulations or quasi three‐dimensional simulations, where flow and transport in aquitards are treated using one‐dimensional analytical and/or numerical approximations. Spatial discretization of three‐dimensional regions is performed using a vertical slicing approach designed to accommodate complex geometry with irregular boundaries, layering, and/or lateral discontinuity. This approach is effectively combined with the use of simple linear elements such as rectangular and triangular prisms, and composite hexahedra and pentahedra made up of tetrahedra. For these elements, computation of element matrices can be performed efficiently using influence coefficient formulas that avoid numerical integration. New transport influence coefficient formulas are presented for rectangular and triangular prism elements. Matrix assembly is performed slice by slice, and the matrix solution is achieved using a slice successive relaxation scheme. This permits a fairly large number of nodal unknowns (of the order of five to ten thousand) to be handled conveniently on small or medium‐size minicomputers. Flexibility of the formulation and matrix handling procedures also allows two‐dimensional and axisymmetric problems to be solved efficiently using single slice representations. Four examples are presented to demonstrate the model verification and utility. These problems represent a fair range of physical conditions. Where possible, simulation results are compared with previously published solutions.