Hantush's well function involves leakage from an overlying layer. The leakage is parameterized by the characteristic length, with the simbol $\lambda$
$$ \lambda = \sqrt{ kD c}$$where $c$ [dimentions time] is the vertical resistance of the overlying aquitard. It can also be expressed as $c = b/k_v$ with $b$ the thickness of the overlying aquitard and $k_v$ its averaged vertical conductivity.
with $\rho = r/\lambda$, Hantush's well function is given by
$$ W_h(u, \rho) = \intop_u^\infty \frac {e^{-y - \frac{ \left( \frac {\rho} 2 \right)^2} {y}}} {y} dy $$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sp
In [17]:
def Wh(u, rho):
'''Hantush's well function (with leakage)
parameters
----------
u : array of floats
rho : float
returns
-------
Wh : array like u with Hantush's well function values
2018-01-12
'''
if np.isscalar(u): # u is a scalar?
u = np.asarray([u]) # if so, turn it into an array
w = np.zeros_like(u) # outcome array, initially all zeros
for i, uu in enumerate(u): # loop ver all values of u array
y = np.logspace(np.log10(uu), 8, 201)
ym = 0.5 * (y[:-1] + y[1:]) # y of midpoints
arg = np.exp(-ym - (rho/2)**2 /ym) / ym # arg at midpoints
dy = np.diff(y) # same as y[1:] - y[:-1]
w[i] = np.sum(arg * dy) # integration
if len(u) == 1: # u was a scalar?
return w[0] # then return first value of array w
else: # otherwise
return w # return the entire array of values
With Wh(u, rho) in place, we can now plot the Hantush type curves, i.e. the graphs of Wh(u, rho) versus 1/u for different values of rho.
In [23]:
u = np.logspace(-6, 1) # suitable range for values of u
Rho = [0.0, 0.01, 0.03, 0.1, 0.3, 1, 3 ] # some values of rho
# notice that for rho=0, W_Hantush = W_Theis
for rho in Rho: # for each rho plot one line
plt.plot(1/u, Wh(u, rho), label='rho={:.2f}'.format(rho))
# embellishment of plot
plt.title('Hantush type curves')
plt.xlabel('1/u')
plt.ylabel('Wh(u)')
plt.xscale('log') # x axis on log scale
plt.yscale('log') # y axis on log scale
plt.legend()
plt.grid()
plt.show() # show it
In [35]:
kD = 500 # m2/d
S = 0.2 # [-]
c = 600 # d
t = 1.2 # d
r = 25. # m
Q = 1200 # m3/d
lamb = np.sqrt(kD * c)
u = r**2 * S / (4 * kD * t)
rho = r/lamb
s = Q/(4 * np.pi * kD) * Wh(u, rho)
print('Drawdown at r={} m and t={} d equals s={:.2f} m'.format(r, t, s[0]))