Embedded Files
Python Program for Analyzing Water Table Aquifer Pumping Test Data
Python Program for Analyzing Water Table Aquifer Pumping Test Data
Introduction to Python Unconfined Aquifer Test Analysis Web Page
Introduction to Python Unconfined Aquifer Test Analysis Web Page
The purpose of this web page is to present a Python program for modern unconfined aquifer pumping test analysis. I developed the Python code by translating and modifying parts of the WTAQ Fortran code developed by Barlow and Moench (1999) and adding (1) an input module and (2) output graphics. The Python program consists of the two modules in Attachment 1 and Attachment 2. The input module is PYWTAQ_INP.py and the main module is PYWTAQ.py. Both modules must be altered when applied to test data from different wells or piezometers. Appropriate input values must be entered in PYWTAQ_INP.py, and the output graph title must be entered in PYWTAQ.py. The program provides a modern efficient way to modify aquifer parameters and instantaneously compare the resulting calculated dimensional time-drawdown curve with the field data curve until the curves match. This web page is technical and assumes knowledge and understanding of water well pumping test analysis and Python programming. The Python program was initially tested by applying it to test data supplied by Barlow and Moench (1999). Then it was applied to field data supplied by Green and others (1991).
Test Data Analysis
Test Data Analysis
Barlow and Moench (1999) supply test data for a simulated well and four piezometers in their Table 5. The simulated parameter values for the aquifer and pumped well are given in their Table3. Using this data in PYWTAQ_INP.py and running PYWTAQ.py provides a comparison of WTAQ test data and PYWTAQ results. Figure 1 shows the comparison for the pumping well. The results are identical, proving that the PYWTAQ code related to the pumping well is correct.
Figure 1. Comparison of results for pumping well test. The WTAQ data are green and the PYWTAQ results are red.
Similarly, the results for a peizometer are compared in Figure 2. The results are close, indicating that the PYWTAQ code related to observation piezometers is acceptable .
Figure 2. Comparison of results for piezometer test.
Application of PYWTAQ to Field Pumping Test Data
Application of PYWTAQ to Field Pumping Test Data
PYWTAQ was applied to analyze actual pumping test data provided by Green and others (1991). The pumping test data is for the pumping well and two observation wells.
Description of Water Table Aquifer Tested
Description of Water Table Aquifer Tested
The water table aquifer that was tested is composed of fine to medium grained silty sand with lenses of fine-grained calcareous sandstone. Its saturated thickness was 98 feet. The aquifer is underlain by nearly impermeable sandstone interbedded with clay.
Pumped Well Description
Pumped Well Description
The pumped well was a 28-inch diameter bore-hole. The screened interval extended from 18 feet below the initial water table to the base of the aquifer. The screened interval was complicated. It consisted of a 12-inch diameter wire-wrapped steel screen inside 16-inch diameter slotted casing. The annulus between the screen and the casing contained a sand pack with a mean diameter of 1/16 inches, and the annulus between the slotted casing and the wall of the hole was packed with 1/4 to 3/8 inch pea gravel.
Observation Well Description
Observation Well Description
Two observation wells were used. They were cased with slotted 2-inch PVC pipe. The hole diameters of the observation wells are not given. The nearest observation well (E-5) was 88 feet from the pumped well and was screened in the upper 58 percent of the saturated aquifer. The other observation well (S-4) was 111 feet from the pumped well and was screened through the entire thickness of the saturated aquifer.
Description of Original Pumping Tests and Analysis
Description of Original Pumping Tests and Analysis
In the original investigation, a step test was performed, and graphical analysis that assumed a value of 2 for the well loss exponent yielded a well loss coefficient of 2.1E-3 ft min2/gal2.
Two constant rate pumping tests were performed. The first constant rate pumping test was conducted at 90 gpm for 42 hours with less than 5 percent discharge variation. Boulton analysis was applied to time-drawdown data from the two monitoring wells. The result was as follows:
Transmissivity: 1250 ft2/day,
Hydraulic conductivity: 12.76 ft/day,
Storativity: 7E-3 (dimensionless),
Specific yield: 0.03 (dimensionless).
The second constant rate pumping test was conducted at 150 gpm, and the results were the same.
Pumping Well Step Test Analysis
Pumping Well Step Test Analysis
PYWTAQ uses a well loss coefficient, which can be obtained from a step test of the pumping well. Green and others (1991) provide step test data. This data was analyzed for the present study. In this analysis "adjusted time" was plotted against "adjusted drawdown" in the manner described on the step test web page. The resulting plot (Figure 3) was subjectively fitted by parallel straight lines for each step. Such straight lines on a semi-logarithmic graph imply that, for the period of the test, the drawdown can be described by an equation of the form s=C1Qlog(C2t), where s is drawdown, t is time since the test began, and C1 and C2 are constants. Consequently, values for the well loss coefficient can be obtained from the vertical separation of the lines. The determination of well loss from such plots is described on the step test web page. The three values obtained from the plot in Figure 1 are 1.3E-3, 1.5E-3, and 1.1E-3 ft min2/gal2 (average 1.3E-3). This well loss is the same order of magnitude as the one obtained in the original graphical analysis cited above ( 2.1E-3 ft min2/gal2). The well loss exponent was assumed to be 2 in the original analysis and in the present analysis.
Figure 3. Step test plot.
PYWTAQ Constant Rate Pumping Test Analysis
PYWTAQ Constant Rate Pumping Test Analysis
The present analysis used PYWTAQ to analyze the data from the 90 gpm pumping test.
PYWTAQ Dimensional-Curve Analysis of Pumped Well Data
PYWTAQ Dimensional-Curve Analysis of Pumped Well Data
Figure 4 shows the dimensional-curve match with the observed drawdown data from the pumped well.
Figure 4. Dimensional-curve match for pumped well. The red curve is the PYWTAQ result and the green curve is the observed data.
The following aquifer coefficients were used to get this fit:
Horizontal hydraulic conductivity (hkr): 8.0 ft/day.
Vertical hydraulic conductivity (hkz): 5.0 ft/day.
Aquifer specific storage (ss): 7.0E-4 (1/ft).
Aquifer specific yield (sy): 0.03 (dimensionless)
The calculated curve for the pumping well assumes (1) instantaneous drainage at the water table, and (2) a well loss coefficient of 4.69E-9 (day2/ft5).
PYWTAQ Dimensional-Curve Analysis of Data for Observation Wells
PYWTAQ Dimensional-Curve Analysis of Data for Observation Wells
Data from two observationwells were analyzed. Figure 5 shows the curve match for well E-5, a partially penetrating well located 88 feet from the pumping well. It penetrated 58 percent of the saturated thickness of the aquifer (98 feet). The diameter of the well was 2 inches, and it was perforated its entire length. The aquifer coefficients used to get this fit are the same as for the pumping well.
Figure 5. Dimensional curve match for observation well E-5.
Figure 6 shows the curve match for observation well S-4. It was a fully penetrating 2-inch well located 111 feet from the pumping well. It was perforated its entire length. The aquifer horizontal hydraulic conductivity and vertical hydraulic conductivity obtained from the curve match are the same as for the pumping well and observation well W-5. However, the specific storage of the aquifer is reduced to 4E-5, and the specific yield is increased to 0.05. Values for specific yield of 0.03 and 0.05 for specific yield are very low. Green and others (1991) suggest that the low specific yield might be partially explained by the large amount of silt and clay in the aquifer. It is noteworty that PYWTAQ does not treat delayed yield at the water table, which might affect the calculated specific yield of the aquifer.
Figure 6. Dimensional curve match for observation well S-4.
Conclusions
Conclusions
PYWTAQ worked well when compared with test data for a pumping well and peizometers. When used to anayze field data for a pumping well and small diameter observation wells, PYWTAQ results are comparable to published results obtained using Boulton analysis.
Attachment 1: PYWTAQ Python Program for Unconfined Aquifer Pumping Test Analysis
Attachment 1: PYWTAQ Python Program for Unconfined Aquifer Pumping Test Analysis
'''
PYWTAQ.py
This is a translation and revision of the Fortran code WTAQ that was developed
by Paul Barlow and Allen Moench, USGS. It generates dimensional curves for
unconfined aquifer constant rate pumping tests. Input is from PYWTAQ_INP. This
version uses day and foot units. Any use of this Python code or derivatives is
the responsibility of the user. Under no circumststances shall the author be
liable for damages that result from such use. This code has not been tested
extensively and may contain errors.
'''
import numpy as np
import math
import sys
from PYWTAQ_INP import *
from scipy.special import k0, k1
import matplotlib.pyplot as plt
#Print to txt file
#f = open('test2.txt', 'w')
#sys.stdout = f
#Function for stehfest coeffiecients:
def stehfest_coefficients(N):
"""
Calculates the Stehfest coefficients for a given number of terms N.
Args:
N: An even integer representing the number of terms.
Returns:
A NumPy array containing the Stehfest coefficients (v_j) for j = 1 to N.
"""
if N % 2 != 0: # % divides N by 2 and returns the remainder.
raise ValueError("N must be an even integer.")
v = np.zeros(N + 1)
for j in range(1, N + 1):
sum_val = 0
for k in range(int((j + 1) / 2), min(j, N // 2) + 1):
numerator = k**(N // 2) * math.factorial(2 * k)
denominator = math.factorial(N // 2 - k) * math.factorial(k) *\
math.factorial(k - 1) * math.factorial(j - k) * math.factorial(2 * k - j)
sum_val += numerator / denominator
v[j] = (-1)**(j + N // 2) * sum_val
return v[1:] # Return coefficients from v[1] to v[N]
#Function for Laplace transform solution for hd.
def ltst3(td, ipws, iows, idpr, ns, weltyp, ntms,
betaw, sigma, rerrnr, v, xln2, expmax,
rd, zd, zd1, zd2, xld, xdd, wd, beta):
"""
Calculates the Laplace transform solution for drawdown.
"""
# --- Initialization ---
hd = 0.0
xn_max = ntms * (2.0**((-np.log10(beta) - 2.0)))
nx = int(round(xn_max))
if nx < 4:
nx = 4
pi = np.pi
# Handle fully penetrating pumped well
if ipws == 1:
xdd = 0.0
xld = 1.0
# Handle fully penetrating observation well
if iows == 1:
zd1 = 0.0
zd2 = 1.0
xp = 0.0
# --- Main Stehfest Loop ---
for i in range(1, ns + 1): #ns is number of Stehfest terms
pp = xln2 * i / td
# --- RHS Calculation ---
rhs = 0.0
rhs = pp / (sigma * betaw)
suma = 0.0
sume = 0.0
#Use Newton-Raphson to calculate epsilon for current value of pp.
nnn = 0
# Initialize eps (root) variable
eps = 0.0
eps0 = 0.0
# --- Finite Sum Loop (Root Finding) ---
while nnn < nx:
nnn += 1
# --- Guess Initial EPS ---
if nnn == 1:
# First root guess
if rhs < 1.0:
eps = np.sqrt(rhs)
else:
eps = np.arctan(rhs) #In radians.
else:
# Subsequent roots are approx PI apart
eps = eps0 + pi
# --- Newton-Raphson Iteration ---
n = 0
while True:
eps0 = eps
n += 1
if n > 100:
print('Program stopped Newton-Raphson iterations exceed 100')
exit()
#Calculate function F and its derifative F prime
a1 = np.sin(eps0)
a2 = np.cos(eps0)
f = eps0 * a1 - rhs * a2
fp = eps0 * a2 + a1 + rhs * a1
# Newton-Raphson update: x_new = x_old - f(x)/f'(x)
eps = eps0 - f / fp
# Check convergence
if abs(eps - eps0) <= rerrnr * abs(eps):
break # Exits the Newton-Raphson loop.
# --- Calculate Terms using Converged EPS ---
qn = np.sqrt(betaw * eps**2 + pp)
if qn > expmax:
qn = expmax
db = np.sin(eps * (1.0 - xdd))
da = np.sin(eps * (1.0 - xld))
if ipws == 1:
da = 0.0
sines = db - da
# Bessel functions K0 and K1
re0 = k0(qn) #Modified Bessel Functions of the second kind order 0
re1 = k1(qn) # Order 1
# --- Calculate 'A' Term (Pumped Well) ---
xnum = re0 * sines**2 / (eps * (xld - xdd))
xden = 0.5 * qn * re1 * (eps + 0.5 * np.sin(2.0 * eps))
a = xnum / xden
suma += a
# --- Calculate 'E' Term (Observation Well) ---
if weltyp > 1:
qnrd = qn * rd
if qnrd > expmax:
qnrd = expmax
re0x = k0(qnrd)
xnum_e = 0.0
if iows == 0:
# Partially penetrating obs well
xnum_e = re0x * sines * ((np.sin(eps * zd2) -
np.sin(eps * zd1)) / (eps * (zd2 - zd1)))
elif iows == 1:
# Fully penetrating obs well
xnum_e = re0x * sines * np.sin(eps) / eps
elif iows == 2:
# Observation piezometer
xnum_e = re0x * sines * np.cos(eps * zd)
term_e = xnum_e / xden
sume += term_e
# End of NNN loop (implied by while)
# --- Aggregate Stehfest Terms ---
denom = (1.0 + wd * pp * suma)
pdl = 0.0
if weltyp == 1:
pdl = suma / (pp * denom)
elif weltyp > 1:
if idpr == 0: #Option for delayed response of obs well.
pdl = sume / (pp * denom)
elif idpr == 1: #Delayed response of observation well.
slugf = 1.0 / (1.0 + wdp * pp)
pdl = slugf * sume / (pp * denom)
# Accumulate XP using Stehfest coefficient V
xp = xp + v[i-1] * pdl
# --- Final Result ---
hd = 2.0 * xp * xln2 / (td * (xld - xdd))
return hd
# Define selected program parameters.
f1 = (ss*rw*rw)/hkr #(days)
f2 = qq/(4.0*math.pi*hkr*bb)
xkd = hkz/hkr
rwd = rw/bb
betaw = xkd * rwd * rwd
expmax = 708.
v = stehfest_coefficients(ns)
xln2 = math.log(2.0)
nts = nox*nlc + 1
deltd = 10.0**(1.0/nox) #(1)
td = tlast*deltd/f1 #(1)
sigma = (ss*bb)/sy
zd=zd1=zd2 = 0.0 # To initiate var
wd = 0.0
sw = 0.0
if weltyp == 1:
zd = 0.
wd = ((rc/rw)**2)/(2.*ss*(zpl-zpd))
if ipws == 0: #Partially penetrating
xdd=zpd/bb #(1)
xld=zpl/bb #(1)
if ipws == 1: #Fully penetrating
xdd = 0.0
xld = 1.0
rd=1.0
rdsq = rd * rd
beta = betaw * rdsq
if weltyp == 2:
rd = r / rw
rdsq = rd * rd
beta = betaw * rdsq
xdd = zpd / bb
xld = zpl / bb
if iows == 0: # Partially penetrating
zd, zd1, zd2 = 0.0, 1.0 - (z2/bb), 1.0 - (z1/bb)
elif iows == 1: # Fully penetrating
zd, zd1, zd2 = 0.0, 0.0, 1.0
if idpr == 0:
wdp = 0.0
elif idpr == 1:
xm = math.sqrt(xkd)
xx = xll / (2.0 * xm * rp)
f_prime = xll / math.log(xx + math.sqrt(1.0 + xx**2))
wdp = ((rp/rw)**2) / (2.0 * ss * f_prime)
if weltyp == 3:
rd = r / rw
rdsq = rd * rd
beta = betaw * rdsq
zd = 1.0 - (zp/bb)
zd1 = 0.0
zd2 = 0.0
xdd = zpd / bb
xld = zpl / bb
ihd, ihdh, init = 1, 1, 1 #Initial values for drawdown calculations.
tdplot = [] #To initiate tdplot.
hdplot = []
for nt in range(1, nts + 1):
td = td/deltd
tdplot.append(td)
hd_result = ltst3(
td=td, # TD: Dimensionless time
ipws=ipws, # IPWS: Partial penetration pumped well flag
iows=iows, # IOWS: Partial penetration obs well flag
idpr=idpr, # IDPR: Delayed response flag
ns=ns, # NS: Number of Stehfest terms
weltyp=weltyp, # weltyp: 0 is pumped, 1 is obs., 3 is peizometer
ntms=30, # NTMS: Max terms for summation. Suggested: 20 or 30
betaw=betaw, # BETAW: Dimensionless beta for well
sigma=sigma, # SIGMA: Specific yield/storage ratio
rerrnr=rerrnr, # RERRNR: Newton-Raphson tolerance
v=v, # V: Stehfest coefficients array
xln2=xln2, # XLN2: ln(2)
expmax=700.0, # EXPMAX: Max exponent limit
rd=rd, # RD: Dimensionless radius
zd=zd, # ZD: Dimensionless vertical coordinate
zd1=zd1, # ZD1: Obs well screen bottom
zd2=zd2, # ZD2: Obs well screen top
# xld=zpl, # XLD: Pumped well screen bottom
# xdd=zpd, # XDD: Pumped well screen top
xld=xld, #Gemini recomendation
xdd=xdd, #Gemini recomendation
wd=wd, # WD: Well bore storage parameter
beta=beta # BETA: Factor for series convergence
)
hdplot.append(hd_result)
tdplotrr = [item*f1*(1/(r*r)) for item in tdplot]
if weltyp == 1:
hdplotdim = [item * f2 + (wc * qq *qq) for item in hdplot]
if weltyp > 1:
hdplotdim = [item * f2 for item in hdplot]
adjt = [item / (1440. * r *r) for item in obst] # (days/(feet*feet))
#Begin plotting.
plt.figure(figsize=(8.0,6.0)) #Adjust figure size (decimal inches) as needed.
plt.xscale('log')
plt.yscale('log')
plt.minorticks_on()
plt.tick_params(which='minor',length=10,width=2,colors="k",direction='inout')
#Get the current axes object
ax = plt.gca() #Sets ax as name of the axis object of the current axes.
#Set the linewidth for each of the spines (top, bottom, left, right)
#Adjust the linewidth to desired thickness
ax.spines['top'].set_linewidth(4.0)
ax.spines['bottom'].set_linewidth(4.0)
ax.spines['left'].set_linewidth(4.0)
ax.spines['right'].set_linewidth(4.0)
#Set the color for each of the spines.
ax.spines['top'].set_color('k') # Color for the top spine
ax.spines['bottom'].set_color('k') # Color for the bottom spine
ax.spines['left'].set_color('k') # Color for the left spine
ax.spines['right'].set_color('k') # Color for the right spine
#Modify this block to get the labels and title wanted.
ax.set_xlabel('TIME DIVIDED BY RADIUS SQUARED (DAYS/(FT*FT))', fontsize = 16)
ax.set_ylabel('DRAWDOWN (FT)', fontsize = 16)
ax.set_title('WATER TABLE AQUIFER TEST (Pine S-4)')
ax.grid()
#Optional: Tick marks thickness.
#plt.tick_params(axis='x', which='both', width=1.5, length=7)
#plt.tick_params(axis='y', which='both', width=1.5, length=7)
plt.plot(tdplotrr,hdplotdim,linewidth=8, color ='red') #adjtd is on x-axis.
plt.plot(adjt,obsdd, linewidth=4, color = 'green') # Plots on same grapth.
plt.savefig('adjhd_vs_adjtd.png')
# Restoring output back to the terminal
#sys.stdout = sys.__stdout__
print('Close the show window to complete execution')
plt.show()
print('End PYWTAQ.py')
Attachment 2: Input Module for PYLAQ
Attachment 2: Input Module for PYLAQ
'''
PYWTAQ_INP.py
Start with estimates of ss, sy, hkr, hkz.
'''
import numpy as np
#Input for Pine Ridge observation well W-5.
bb = 98.0 #Initial saturated thickness of the aquifer. (ft)
hkr = 8.0 #Horizontal hydraulic conductivity of the aquifer. (ft/day)
hkz = 5.0 #Vertical hydraulic conductivity of the aquifer. (ft/day)
ss = 4.0e-5 #Specific storage of the aquifer. (1/ft)
sy = 0.05 #Specific yield of the aquifer. (1)
tlast = 1.75 #Largest value of time. (days)
nlc = 3 #Number of log cycles on the time scale for which drawdown (dd) calculated. (1)
nox = 3 #Number of equally spaced times per log cycle for which dd calculated. (1)
rerrnr = 1.0e-10 #Relative error for Newton-Raphson iteration. 1.0e-10 suggested. (1)
ns = 8 # Number of Stehfest terms. Even number. Eight usually sufficient. (1)
weltyp = 2 #Type of well. Pumped=1. Observation=2. Piezometer=3.
ipws = 1 #Type of pumped well. 0 = partially penetrating. 1 = fully penetrating. (1).
qq = 17325 #Pumping rate. (ft3/day)
rw = 0.5 #Radius of pumped well screen. (ft)
rc = 1.17 #Inside radius of pumped well where water level is changing. (ft)
zpd = 18.0 #Depth from initial water table to top of screen(s) of pumped well. (ft)
zpl = 98.0 #Depth from initial water table to bottom of screen(s) of pumped well.
wc = 4.6897e-9 #Well loss coefficient. (day^2/ft^5)
iows = 1 #Type of observation well. Partially penetrating = 0. Fully penetrating = 1.
#Piezometer = 2 (1)
idpr = 1 #Option for delayed response of observation well. Does not apply to a piezometer.
#Set to 0 for a peizometer. No delayed response = 0. Delayed response = 1.
r = 111.0 #Radial distance from axis of pumped well to observation location.
#Use 1.0 for the pumped well (ft)
z1 = 0. # Depth below initial water table to top of screened interval of obs well. Piezometer = 0.(ft)
z2 = 98.0 # Depth below top of aquifer or initial water table to bottom of screened interval.
#Set to bb for fully pentrating, Piezometer = 0. (ft)
zp = 0.0 #Depth below top of aquifer or initial water table to center of piezometer. Well = 0. (ft)
rp = 0.0833 #Inside radius of the observation well in the interval where water level
# is changing. 0.0 if idpr = 0.(ft)
xll = 98.0 #Length of screened interval of observation well. 0 if idpr = 0. (ft)
# Read TIMEOB and XMEASoB from text file. Program should plot calculated and observed drawdowns.
column1, column2 = np.loadtxt('Pine_S-4.txt', unpack=True)
obst = column1.tolist() #Converts the numpy array to a list.
obsdd = column2.tolist()
print('End PYWTAQ_INP.py')
Attachment 3: Drawdown Text File for PYWTAQ_INP
Attachment 3: Drawdown Text File for PYWTAQ_INP
7 0.01
9 0.02
10 0.03
12 0.04
15 0.06
20 0.1
30 0.18
40 0.29
50 0.38
60 0.49
70 0.55
80 0.59
90 0.62
100 0.65
120 0.7
150 0.76
200 0.84
302 0.97
400 1.13
504 1.28
605 1.42
705 1.54
806 1.6
905 1.7
1008 1.8
1215 1.97
1382 2.07
1496 2.14
1729 2.29
1886 2.37
2020 2.47
2255 2.62
2497 2.73
Bibliography for PYWTAQ Leaky Aquifer Pumping Test
Bibliography for PYWTAQ Leaky Aquifer Pumping Test
Barlow, Paul M. and Allen F. Moench (1999): WTAQ--A Computer Program for Calculating Drawdowns and Estimating Hydraulic Properties for Confined and Water-Table Aquifers; U.S. Geolgical Survey Water-Resources Investigations Report 99-4225. 🔗
Green, E. A. and others (1991): Aquifer Tests and Water-Quality Analyses of the Arikaree Formation Near Pine Ridge, South Dakota; U.S. Geological Survey Water-Resources Investigations Report 91-4005. 🔗
Bibliography for PYWTAQ Leaky Aquifer Pumping Test
Bibliography for PYWTAQ Leaky Aquifer Pumping Test
Barlow, Paul M. and Allen F. Moench (1999): WTAQ--A Computer Program for Calculating Drawdowns and Estimating Hydraulic Properties for Confined and Water-Table Aquifers; U.S. Geolgical Survey Water-Resources Investigations Report 99-4225. 🔗
Green, E. A. and others (1991): Aquifer Tests and Water-Quality Analyses of the Arikaree Formation Near Pine Ridge, South Dakota; U.S. Geological Survey Water-Resources Investigations Report 91-4005. 🔗
Report abuse