Finite difference and finite element methods are common numerical techniques for modeling groundwater flow problems. These methods differ from those discussed thus far in that the hydraulic head is approximated by discretization throughout the flow domain. As a result, not only are the boundary conditions approximated, but also the differential equation itself. Characteristic for both methods is that the flow domain is bounded. This is true also for applications to regional flow in contrast to the analytic element method, where the aquifer system is modeled as being infinite in extent; boundary conditions are applied only along internal boundaries that are physically present in the aquifer system.

Advantages of both the finite difference and finite element methods are that the hydraulic conductivity can be easily varied throughout the aquifer system, the formulations are well suited for modeling transient flow, and they are comparatively straightforward. Both methods are discussed in detail in the literature on groundwater mechanics, e.g., Verruijt [1982], Wang and Anderson [1982], and Bear and Verruijt [1987].

A number of computer implementations of finite difference and finite element techniques exist. The most popular finite difference program is MODFLOW, e.g., Mehl and Hill [2001], [2003], [2005], and MODFLOW-2005. A popular finite element program is FEFLOW; manuals for this program are available online.

There is an extensive literature on groundwater flow modeling using numerical methods.We limit the treatment in this text to a discussion of the basic principles of such models; it is important, when applying numerical software, to understand the basic principles involved, as well as the strengths and weaknesses of the methods; we focus on these aspects.

Finite Difference Methods

Finite difference methods are based on a discretization of the flow domain into a mesh that is usually rectangular, where the potential or head is computed at the gridpoints by solving the differential equation in finite difference form throughout the mesh. There exist a variety of numerical techniques for solving the resulting system of linear equations, and the method may be applied with comparative ease to problems where the hydraulic conductivity varies from node to node. Our purpose is to explain the basic principles involved in the method for a square mesh and constant hydraulic conductivity.