PYTHON FINITE-DIFFERENCE GROUNDWATER MODEL

A Compact Computer Program for a Simple Numerical Groundwater Model

By Darrel Dunn, Ph.D., Hydrogeologist

Purpose and Scope of Python Groundwater Model Webpage

The purpose of this webpage is to present python code for a simple steady-state finite-difference groundwater model. The presentation is moderately technical, requiring basic knowledge of the Python language, algebra, and groundwater hydrology. Python has become widely used by scientists and engineers because it is concise, readable and powerful. Many useful libraries, modules and functions are available. Libraries useful in geology and hydrology include NumPy and Matplotlib. Both are used in finite-difference groundwater modeling. The variables associated with the finite-difference grid can be put in NumPy arrays and manipulated concisely in a variety of ways using NumPy functions. The Matplotlib Pyplot module is useful in plotting the results of the modeling using very few lines of code. In the past, Python and other interpreted languages were impractical for finite-difference modeling because they execute relatively slowly compared to compiled languages like Fortran. In the present, even desktop computers can be designed for great speed Therefore, Python can be used for compatibly sized finite-difference models. Furthermore, Python can interface with compiled code for performance enhancement.

This webpage presents Python code for a very simple finite-difference model, just to illustrate how Python can be used. Such a simple model can be expanded and modified to satisfy the objectives of a particular modeling task. Rather than use complex open-source or proprietary compiled code designed for universal modeling application, code can be written in Python optimally designed for a particular task.

Finite-Difference Equation Using NumPy Arrays

Finite-difference groundwater modeling computer programs such as FDPY.py shown in Attachment 2 are based on an array of cells. Mass-balance for each cell requires that the flow into each cell must equal the flow out in a steady-state system. Consider the cell illustrated in Figure 1. It is for a 1-layer system and the lateral flows out across cell faces are QR, QL, QT, and QB. Another flow out is depicted as Qout, which might be groundwater extraction or recharge. Flow in is negative flow out. These must balance:

QR + QL + QT +QB = Qout. (Equation 1)

Diagram of a cell in a Python finite-difference groundwater model.

Figure 1. Cell [1,1] in NumPy array.

In the following mass-balance equations:

The variable dim is the dimensions of the cell in both the x and y directions in feet (cells are square and equal),

The variable height is height of the cells (all heights are equal).

The variables r and c are NumPy row and column numbers.

KX[r-1,c] is the hydraulic conductivity at the left face of cell[r,c] (gallons per day per square feet).

KX[r-1,c-1] is the hydraulic conductivity at the right face of cell[r,c].

KY[r-1,c-1] is the hydraulic conductivity at the top face of t cell[r,c] (NumPy terminology).

KY[r,c-1] is the hydraulic conductivity at the bottom face of cell[r,c].

The variable qout[r-1,c-1] is flow out of cell[r,c] in gallons per day (gpd).

All of the boundary cells are treated as constant head cells. Consequently, the shapes of KX and KY are (rnum-1,cnum-1), where rnum and cnum are the number of rows and columns in the array of cells. The shape of Qout is (rnum-2,cnum-2).

Applying Darcy's law to Equation 1 for each interior cell:

(KX[r-1,c] * height * {head[r,c+1] - head[r,c]}) + (KX[r-1,c-1] * height * {head[r,c-1] - head[r,c]})

+(KY[r-1,c-1] * height * {head[r-1,c] - head[r,c]}) + (KY[r,c-1] * height * {head[r+1,c] - head[r,c]}) = qout[r-1,c-1]

(Equation 2)

Rearranging this equation so that head[r,c] is on the left side:

head[r,c] = (-qout/height) + (KX[r-1,c] * head[r,c+1]) +

(KX[[r-1,c-1] * head[r,c-1]) +

(KY[r-1,c-1] * head[r-1,c]) +

(KY[r,c-1] * head[r+1,c])/

(KX[r-1,c]) + (KX[r-1,c-1]) + (KY[r-1,c-1]) + (KY[r,c-1] ))

(Equation 3)

Applied to all interior cells in a model, this equation is a set of linear equations in which the number of unknowns (head[r,c]) is equal to the number of equations. Consequently, the system of simultaneous equations can be solved. However, most groundwater models have a large number of cells, so the number of equations is large and a direct solution is cumbersome. The solution can be simplified by using the Gauss-Seidel iteration method. In this method, head[r,c] is successively calculated for each interior cell using previously calculated values. The calculation iterates through all active cells until the heads are only slightly changing within a prescribed difference limit. These final heads are the solution to the system of equations and are the heads that satisfy the mass balance requirement of the model. Gauss-Seidel iteration is guaranteed to converge on a solution if (1) head[r,c] in Equation 3 is treated as the diagonal term in the set of equations and (2) the absolute value of the coefficient of this diagonal is greater than or equal to the sum of the absolute values of the off-diagonal coefficients in each equation, plus (3) the absolute value of the coefficient of one diagonal term is greater than the sum of the absolute values of the off-diagonal coefficients (McCracken and Dorn, 1964). Only the first two requirements are satisfied in Equation 3. Consequently convergence is not guaranteed, nevertheless convergence is usual.

Simple Example of a Python Finite-Difference Groundwater Model

Attachment 2 is a python program (FDPY.py) for a simple one-layer model as described above. It receives its input from FDPY_in. The example of FDPY_in shown in Attachment 1 is for a very simple model that just serves as an illustration of how Python finite-difference modeling works. The model has the following attributes:

Five rows and five columns of cells.
Initial hydraulic heads are assigned to each cell (feet): head.shape is (5,5).
All outer boundary cells are held constant at the initial value.
An array is used to identify interior cells with constant heads with a 1 and variable heads as zero: typ.shape is (3,3).
An array for groundwater extraction (positive) and recharge (negative) for interior cells (gpd); qout:shape is (3,3)
Arrays of hydraulic conductivity at cell faces of interior cells (gpd/square foot): KX.shape is (4,4) and KY.shape is (4,4).
The cell lateral dimensions (dim) are all 100 feet (square).
The cell heights are all 10 feet.
The maximum head difference between iterations for convergence was set at 0.0001 ft.

These values are set in FDPY_in, as shown in Attachment 1. If the model were scaled up, the necessary changes would be made in FDPY_in. FDPY_in contains a relaxation coefficient (r1), which is set at 1 for the Gauss-Seidel method. This coefficient can be set to a value between 1 and 2. For values greater than one, the method is called successive over-relaxation (SOR). In some cases SOR may decrease the number of iterations for convergence and reduce computation time. However, SOR has a greater potential for divergence that depends on the choice of r1.

In this simple example, an extraction rate of 10,000 gpd was set in the center cell (qout[1,1]). Using the Gauss-Seidel method, the computation converged at 15 iterations. Total execution time was 0.073 seconds on a modern desktop computer with a Linux operating system. The resulting steady-state heads are shown in Figure 2, which was produced by head_plot.py. The code for head_plot.py is shown in Attachment 3. The mass-balance was analyzed using FDPY_Q.py, which is shown in Attachment 4. FDPY_Q.py calculates net inflows to the individual cells. The maximum absolute value of net inflow to any cell was 0.28 gpd. The sum of all net inflows was -0.79 gpd. A map of net inflows is shown in Figure 3.

Conclusion Regarding Python Finite-Difference Groundwater Modeling

These results demonstrate that much larger and more complex Python finite-difference groundwater models designed to address specific problems are feasible.

Hydraulic head in Python groundwater model.

Figure 2. Hydraulic head (feet) in a simple Python finite difference steady-state groundwater model with constant boundaries of 1000 feet, and extraction of 10,000 gpd from the central cell.

Net inflow to cells in a Python finite-difference groundwater model.

Figure 3. Net inflow (gpd) to cells in a simple Python finite difference steady-state groundwater model with constant boundaries of 1000 feet, and extraction of 10,000 gpd from the central cell.