Python Finite-Difference Aquifer Pumping Test Analysis

Karstic Aquifer Pumping Test Analysis Using a Python Finite-Difference Model

By Darrel Dunn, Ph.D., Hydrogeologist

Introduction to Python Finite-Difference Aquifer Pumping Test Webpage

The purpose of this webpage is to demonstrate the use of Python finite-difference steady-state groundwater modeling code to analyse steady-state pumping test data.  (The pumping test is long enough for drawdown to be treated as steady.)  This is a technical webpage and assumes the reader knows pumping test procedure and the Python programming language.   Python finite-difference code was introduced in the webpage titled "Python Finite-Difference Groundwater Model."  The code uses Gauss-Seidel itertion to solve a huge array of simultaneous equations.  A relaxation factor may be applied to speed convergence (PSOR).  The basic code presented on that page is modified for the present model to make it run faster.  This modification uses the JIT compiler from the Numba Library to speed the iteration process.  The use of Numba JIT is illustrated in the code presented in Attachment 1.  The pumping test used in this demonstration is the same one analyzed by using DP_LAQ and described in the webpage titled "Fractured Aquifer Testing."  The original pumping test and analyis was conducted by the U.S. Geological Survey (Crouch, 1994).  The original analysis used the Jacob method (Freeze and Cherry, 1979, p.347).  Crouch opined that a numerical simulation of the system might produce improved estimates of aquifer properies and the effects of pumping.  This webpage describes such a numerical simulation of the system.

Here is an overview of the questions answered explicity or implicitly in this webpage:

Can Python be used for analysis of aquifer pumping tests using finite-difference groundwater modeling?

Answer: Yes.  This webpage provides an example of a steady state analysis.

What is the best pumping test analysis method for a karstic aquifer?

Answer:  A finite-difference groundwater model that can simulate inhomogeneous and anisotropic aquifers with highly conductive groundwater flow paths.

Can a Python groundwater model with a very large number of cells execute quickly?

Answer: Yes.  A Python compiler module can be applied to the function that performs the iterative solution.

Does Python offer efficient management of input files during model calibration?

Answer: Yes.  The model grid can be represented by a NumPy array, and hydraulic conductivities can be successively modified by grid slicing during model calibration.  Descriptions of the input variables can be placed as comments next to the variables in the input file.

Can a finite-difference groundwater model produced by aquifer test analysis be used to estimate drawdown throughout the groundwater system?

Answer:  Yes.  The drawdown in any model cell is available and Python modules can produce maps customized to fit project objectives.

Can a finite-difference groundwater model produced by aquifer test analysis be used to determine the best locations for future water wells.

Answer: No.  It only produces effective hydraulic conductivity throughout most of the area.  It does not identify high conductivity locations except near the pumping well of the aquifer test.  However, in this case it suggests locating test wells on topgraphic lineaments.

Does effective Python finite-difference analysis of pumping test data require special knowledge?

Answer: Yes.  Effective analysis requires knowledge of Python coding, computer terminal command language, aquifer pumping test theory, and relevant geologic subjects such as stratigraphy.

Figure 1.  Location of the pumping well (ZS-3) and monitoring points in the central part of the model grid.

Calibration of the Finite-Difference Model

The model was calibrated by successive trials with different hydraulic conductivities (kx, ky) at cell faces.  Trials were continued until simulated drawdowns were accepted to be close enough to measured steady-state drawdowns.  Groundwater flow in the model area is through a succession of aquifer materials of widely different hydraulic conductivities.  However, the exact locations of the high conductivites (caused by solution of limestone) and locations of lower conductivies of less altered rock is unknown.  Consequently, equivalent hydraulic conductivity was used for most of the model area.  Equivalent hydraulic conductivity is discussed by Freeze and Cherry (1979, p. 33).   The resulting ky values are shown in Figure 2.  Figure 3 zooms to the area containing the observation points.  Note the high y-direction hydaulic conductivies following the lineament direction though the pumping well and including ZS-1 and ZS-13.  Figure 4 shows the hydraulic conductivity in the x-direction, which is normal to the lineament direction.  No high hydraulic conductivities were assigned to x-direction cell faces.  They were left at equivalent hydraulic conductivity values.  To see the actual hydraulic conductivites at all cell faces, look at the array slices for kx and ky in the input file.  The degree of calibration that was accepted for each monitoring point is shown in Table 1.

Figure 2.  Hydraulic conductivity in the y-direction in the entire model.  This map was produced by kx_ky_map.py, which is listed in Attachment 3.

Figure 3.  Y-direction hydraulic conductivity in the area of the monitoring points.

Figure 4.  Hydraulic condictivity in the x-direction in the entire model.  This map was produced by kx_ky_map.py, which is listed in Attachment 3.

Table 1.  Measured and calculated drawdown at the monitoring points.

Figure 5 is a map of the hydraulic head in the model area.  The initial head throughout the area was 6312 feet, which remains the head at the constant head bundaries after calibration.

Figure 5.  Hydraulic head in the model area.  This map was produced by head_plot.py, which is listed in Attachment 4.

Convergence of the model required 50216 iterations through the 200,794 cells.  This large number of iterations was required because it was determined through successive trials that the maximum head difference between iterations (mhd) had to be reduced to 0.0000005 feet before the solution stabilized.   Even so, the total time to execute a model run was only 1 minute and 29 seconds using a relaxation coefficient (r1) of 1.2.  The model was run on a desktop computer with an AMD Ryzen 5 8500G processor and a Linux POP!_OS operating system.

Conclusion to Python Finite-Difference Aquifer Pumping Test WebPage

This webpage demonstrates that Python finite-difference groundwater modeling code can be used to analyse pumping test data.  The demonstration applies to pumping tests where the drawdown at the end of the test may be treated as steady-state.  The number of cells in the model can be very large.  A two-dimensional model was used for this demonstration.  The use of Python programs allows for pumping test analysis and presentation of results that may be customized to fit the needs of a particular project by modifying the code.

References

Crouch, T. M. (1991): Evaluation of the Bidahochi and San Andres-Glorieta Aquifers on Parts of the Zuni Indian Reservation, McKinley and Cibola Counties, New Mexico; U.S. Geological Survey Water-Resources Investigations Report 89-4192. 🔗

Crouch, T. M. (1994): Hydrogeology and Aquifer Test of the San Andres-Glorieta Aquifer on the  Southwest Part of the Zuni Indian Reservation, Cibola County, New Mexico; U.S. Geological Survey Water-Resources Investigations Report 94-4033.  🔗

Freeze, R. A. and J. A. Cherry (1979): Groundwater; Prentice-Hall.

Attachment 1:  FDPYJIT.py Python Program for Finite-Difference Groundwater Modeling with JIT Compiler

'''

FDPYJIT.py 

Run with python3. 

Uses gauss-sidel iteration or SOR to simulate a groundwater flow system.   

Author is DarrelDunn.

FDPY_in.py is a script that generates input values used in this program. 

The cells must be square (eqaal dimensions in the x and y directions).

Befor running this program you must have numba installed.  If it is not

installed, run the following commands from within the folder :

sudo apt update

sudo apt install python3-venv python3-pip

python3 -m venv my_env

source my_env/bin/activate

pip install numba

After this, when you open the terminal to run FDPYJIT, run

source my_env/bin/activate

When you are finished working with FDPYJIT and want to leave the isolated 

environment, type deactivate.

This code has not been tested extensively and may contain errors.

Copyright (C) 2026 Darrel Eugene Dunn.  All rights reserved.

For licensing inquiries, or contracting application of this program, contact

ddunn@dunnhydrogeo.com.

'''

import numpy as np

import FDPY_in # Runs FDPY_in and imports its variables.

from numba import njit

@njit #No-PYthon Just-in-Time (JIT) compiler.

def head_iter(nprnum,npcnum,typ,head,KX,KY,qout,height,relax,limit,mx):

    it = 0

    while it < limit:

        maxdiff = 0.

        it = it + 1       

        rcount = 1

        while rcount <= nprnum - 1:  # rcount is row.

            ccount = 1

            while ccount <= npcnum - 1: # ccount is column. 

                if typ[rcount-1, ccount-1] < 1:

                    temp_head=(((head[rcount,ccount+1] * KX[rcount-1,ccount]) +

                              (head[rcount,ccount-1] * KX[rcount-1,ccount-1]) +

                              (head[rcount-1,ccount] * KY[rcount-1,ccount-1]) +

                              (head[rcount+1,ccount] * KY[rcount,ccount-1]) -

                              (qout[rcount-1, ccount-1] / height)) /   

                              (KX[rcount-1, ccount] + KX[rcount-1, ccount-1] +

                               KY[rcount-1, ccount-1] + KY[rcount, ccount-1]))                

                    diff = abs(head[rcount, ccount] - temp_head)

                    if diff > maxdiff: 

                        maxdiff = diff                    

                    # Successive Over-Relaxation

                    head[rcount, ccount] = (relax * temp_head + (1 - relax) *   

                                            head[rcount, ccount])                

                ccount += 1

            rcount += 1           

        print('maxdiff:', maxdiff)

        if maxdiff < mx:

            print('maxdiff < mx reached.', maxdiff, mx)

            print('iterations:', it)

            return head

    return head # Fallback return if limit is hit

# Main Program

cnum = FDPY_in.cnum

npcnum = cnum - 1 # changes to numpy indexes.

rnum = FDPY_in.rnum

nprnum = rnum - 1

relax = FDPY_in.r1

mx = FDPY_in.mhd

limit = FDPY_in.limit

typ = FDPY_in.typ

height = FDPY_in.height

head = FDPY_in.head

KX = FDPY_in.kx 

KY = FDPY_in.ky

qout = FDPY_in.qout

# BEFORE NUMBA RUNS:  Safely verify/cast the data type to floats so Numba 

#doesn't crash on an object array.

head = np.asarray(head, dtype=np.float64)

# Execute JIT accelerated loop.

head =head_iter(nprnum,npcnum,typ,head,KX,KY,qout,height,relax,limit,mx)

np.save('head_binary.npy', head) # Save clean binary file of head array.

print('END FDPY.py')

Attachment2.  FDPY_in_good.  Input File for Calibrated Run of PYFTJIT.py

'''

This program generates input values and arrays for FDPY.py and FDPYJIt.

rnum is the actual number of rows in the model.

cnum is the actual number of columns in the model.

r1 is relaxation coefficient (>=1 and <2)

mhd is maximum difference in head values for convergence.

typ is an array with 1 for constant head values and 0 for values to be calculated.

limit is the number of iterations allowed.

kx and ky are hydraulic conductivity in x and y directions (ft/day)

qout is flow out of each active cell(ft3/day).  Negative qout is inflow.

dim is the lateral dimension of the cells.  It is the same for all cells. (ft)

height is cell height.  It is the same for all cells.  (ft)

'''

import numpy as np

rnum = 450 #Actual number of rows. 

cnum = 450 #Actual number of columns  

r1 = 1.2

mhd = 0.0000005

limit = 100000

typ = np.zeros([rnum-2,cnum-2], dtype=np.int32) 

head = np.full((rnum,cnum), 6312.) #Fig.2 94-033

qout = np.full([rnum-2,cnum-2], 0.)

qout[223,223] = 496650. #qout is ft3/day. The shape is (rnum-2,cnum-2) 

kx = np.full([rnum-1,cnum-1], 1500.)

ky = np.full([rnum-1,cnum-1], 1500.)

ky[0:224,0:450]=20.0 #Model above pumping well.

kx[0:223,0:450]=20.0 #Model above pumping well.

ky[0:223, 223]=5000. #Solution channel from ZS-1 to upper boundary.

ky[225:305, 223] = 1000. #Between ZS-3 and ZS-13+ for solution channels.

dim = 200.

height = 400.

print('END FDPY_in.py')

Attachment 3.  kx_ky_maps.py.  Prints maps of kx and ky 

'''

kx_ky_maps.py

Author: Darrel Dunn

Creates maps of kx and ky, the hydraulic conductivities that are input to PYFDJIT.py to develop a finite-difference groundwater model.

'''

import sys

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.colors as colors

# Import variables directly from FDPY_in.py, which must be in same directory.  

import FDPY_in as pin

def plot_conductivity_maps(kx, ky, dim, save_prefix="conductivity"):

    """

    Generates and saves PNG color maps for kx and ky arrays with colorbars.

    Parameters:

    - kx, ky: 2D numpy arrays of hydraulic conductivity.

    - dim: Lateral dimension of the cells (in feet) for axis scaling.

    - save_prefix: Prefix for the output PNG filenames.

    """

    # Get shape from the imported arrays

    ny, nx = kx.shape

    extent = [0, nx * dim, 0, ny * dim]#Length of axes in feet.

    # Determine global min and max for a shared logarithmic colorbar scale

    vmin = min(kx.min(), ky.min())

    vmax = max(kx.max(), ky.max())

    # Use LogNorm to handle the wide dynamic range (e.g., 10 to 100,000)

    norm = colors.LogNorm(vmin=max(vmin, 1e-3), vmax=vmax)

    cmap = 'viridis'  # Use viridis colorscheme.

    # Pack data in dictionary datasets for clean iterative plotting 

    datasets = {

        'kx': (kx, 'X-direction Hydraulic Conductivity ($K_x$)'),

        'ky': (ky, 'Y-direction Hydraulic Conductivity ($K_y$)')

    }

    #Loop through kx, ky:

    for key, (data, title) in datasets.items():#key captures kx and ky.

        plt.figure(figsize=(8, 7), dpi=300)  # High DPI for crisp PNG graphics

        # Origin 'upper' causes matrix indexing to start at top-left.]

        # extent = extent causes th axes to be labelled from the lower-left.

        im = plt.imshow(data, cmap=cmap, norm=norm, extent=extent, origin='upper')

        # Labels and styling.

        #'title' comes from the datasets dictionary. 

        plt.title(title, fontsize=14, pad=15)

        plt.xlabel('Distance (ft)', fontsize=11)

        plt.ylabel('Distance (ft)', fontsize=11)

        # Colorbar setup

        #fraction is width of color bar as a fraction of the main plot's width.

        #pad defubes the width of the space between the main plot and the bar.

        cbar = plt.colorbar(im, fraction=0.046, pad=0.04)

        cbar.set_label('Conductivity (ft/day)', fontsize=11)

        #Clean up spacing in the figure to make it look nice. 

        plt.tight_layout()

        # Save output

        filename = f"{save_prefix}_{key}.png"

        plt.savefig(filename, bbox_inches='tight')

        plt.close()

        print(f"Successfully generated: {filename}")

# Execution block of script

print("Importing grid parameters from FDPY_in.py...")

# Safely pull the properties generated by FDPY_in.py

# None is fallback if kx does not exist.

kx_data = getattr(pin, 'kx', None)

ky_data = getattr(pin, 'ky', None)

dim_val = getattr(pin, 'dim', 200.)

# Generate the color maps

plot_conductivity_maps(kx_data, ky_data, dim=dim_val)

print('End kx_ky_maps')    

Attachment 4.  head_plot.py.  Prints Maps Hydrauic Head and Drawdown produced by PYFDJIT.py.

'''

head_plot.py

Author: Darrel Dunn

Plots grid of head values from head_value.npy, which is output from PYFDJIT.py.

'''

import numpy as np

from pylab import *

import FDPY_in

cnum = FDPY_in.cnum

npcnum = cnum -1 #changes to numpy indexes.

rnum = FDPY_in.rnum

nprnum = rnum-1

head_array = np.load('head_binary.npy', allow_pickle=True)

#print(head_array)

imshow(head_array, interpolation='none')

plt.colorbar()

plt.savefig('head_array.png')

plt.clf() #Clears the figure so the following one does not plot on top of it.

initial_head = np.full((rnum,cnum), 6312.)

#print('shape initial_head, head_array',initial_head.shape, head_array.shape)

drawdown = initial_head - head_array

#print('drawdown',drawdown)

imshow(drawdown, interpolation='none')

plt.colorbar()

plt.savefig('drawdown.png')

print('Observation Point DD vs Calculated DD')

print('ZS-13,  1.13',drawdown[226,224])

print('ZS-1,   0.90',drawdown[223,224])

print('ZS-101, 0.46',drawdown[251,203])

print('ZS-100, 0.31',drawdown[282,167])

print('RS,     0.15',drawdown[302,121])

print('ZS-3 pumping',drawdown[224,224])

print('END')