In the previous section the Theis well function was introduced. The function, which is in fact the function known as exponential integral by mathematicians, proved available in the standard library of Ptyhon module scipy.special. We modified it a little to make it match the Theis well function exactly and gave it the name "W" like it has in groundwater hydrology books. Then we used it in some examples.
In this chapter we will investigate the Theis well funchtion character a more accurately.
Instead of looking for the function in the available library we could have computed the function ourselfs, for instance by numerical integration.
$$ W(u) = \intop_u^{-\infty} \frac {e^{-y}} y dy \approx \sum_0^N \frac {e^{-y_i}} {y_i} \Delta y_i $$where $y_0 = u_0$ and $N$ has a a sufficiently large value.
In [48]:
import scipy.special as sp
In [49]:
import numpy as np
from scipy.special import expi
def W(u): return -expi(-u)
def W1(u):
"""Returns Theis' well function axpproximation by numerical intergration
Works only for scalar u
"""
if not np.isscalar(u):
raise ValueError("","u must be a scalar")
LOG10INF = 2 # sufficient as exp(-100) is in the order of 1e-50
y = np.logspace(np.log10(u), LOG10INF, 1000) # we use thousand intermediate values
ym = 0.5 * (y[:-1] + y[1:])
Dy = np.diff(y)
w = np.sum( np.exp(-ym) / ym * Dy )
return w
Try it out
In [50]:
U = 4 * 10** -np.arange(11.) # generates values 4, 4e-1, 4e-2 .. 4e-10
print("{:>10s} {:>10s} {:>10s}".format('u ', 'W(u)','W1(u) '))
for u in U:
print("{0:10.1e} {1:10.4e} {2:10.4e}".format(u, W(u), W1(u)))
Is seems that our numerical integration is a fair approximation to four significant digits, but not better, even when computed with 1000 steps as we did. So it is relatively easy to create one's own numerically computed value of an analytical expression like the exponential integral
The theis well function can be expressed also as a power series. This expression has certain advanages as it gives insight into the behavior of its character and allows important simplifications and deductions.
$$ W(u) = -0.5773 - \ln(u) + u - \frac {u^2} {2 . 2!} + \frac {u^3} {3 . 3!} - \frac {u^4} {4 . 4!} + ... $$This series too can be readily numerially comptuted by first defining a function for it. The sum will be computed in a loop. To prevent having to compute faculties, it is easiest to compute each successive term from the previous one.
So to get from term m to term n+1:
$$ \frac {u^{n+1}} {(n+1) . (n+1)!} = \frac {u^n} { n . n!} \times \frac {u \, n} {(n+1)^2} $$This series is implemented below.
In [51]:
def W2(u):
"""Returns Theis well function computed as a power series"""
tol = 1e-5
w = -0.5772 -np.log(u) + u
a = u
for n in range(1, 100):
a = -a * u * n / (n+1)**2 # new term (next term)
w += a
if np.all(a) < tol:
return w
Compare the three methods of computing the well function.
In [52]:
U = 4.0 * 10** -np.arange(11.) # generates values 4, 4e-1, 4e-2 .. 4e-10
print("{:>10s} {:>10s} {:>10s} {:>10s}".format('u ', 'W(u) ','W1(u) ', 'W2(u) '))
for u in U:
print("{0:10.1e} {1:10.4e} {2:10.4e} {2:10.4e}".format(u, W(u), W1(u), W2(u)))
We see that all three methods yiedld the same results.
Next we show the well function as it shown in groundwater hydrology books.
In [53]:
u = np.logspace(-7, 1, 71)
import matplotlib.pylab as plt
fig1= plt.figure()
ax1 = fig1.add_subplot(111)
ax1.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve versus u', yscale='log', xscale='log')
ax1.grid(True)
ax1.plot(u, W(u), 'b', label='-expi(-u)')
#ax1.plot(u, W1(u), 'rx', label='integal') # works only for scalars
ax1.plot(u, W2(u), 'g+', label='power series')
ax1.legend(loc='best')
plt.show()
The curve W(u) versus u runs counter intuitively which and is, therefore, confusing. Therefore, it generally presented as W(u) versus 1/u instead as shown below
In [54]:
fig2 = plt.figure()
ax2 = fig2.add_subplot(111)
ax2.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve versus 1/u', yscale='log', xscale='log')
ax2.grid(True)
ax2.plot(1/u, W(u))
plt.show()
Now W(u) resembles the actual drawdown, which increases with time.
The reason that this is so, becomes clear from the fact that
$$ u = \frac {r^2 S} {4 kD t} $$and that
$$ \frac 1 u = \frac {4 kDt} {r^2 S} = \frac {4 kD} S \frac t {r^2} $$which shows that $\frac 1 u$ increases with time, so that the values of $\frac 1 u$ on the $\frac 1 u$ axis are propotional with time and so the drawdown, i.e., the well function $W(u)$ increases with $\frac 1 u$, which is less confusing.
The graph of $W(u)$ versus $\frac 1 u$ is called the Theis type curve. It's vertical axis is proportional to the drawdown and its horizontal axis proportional to time.
The same curve is shown below but now on linear vertical scale and a logarithmic horizontal scale. The vertical scale was reversed (see values on y-axis) to obtain a curve that illustrates the decline of groundwater head with time caused by the extraction. This way of presending is probably least confusing when reading the curve.
In [55]:
fig2 = plt.figure()
ax2 = fig2.add_subplot(111)
ax2.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve versus 1/u', yscale='linear', xscale='log')
ax2.grid(True)
ax2.plot(1/u, W(u))
ax2.invert_yaxis()
plt.show()
We see that after some time, the drawdown is linear when only the time-axis is logarithmic. This suggests that a logarithmic approximation of time-drawdown curve is accurate after some time.
That this is indeed the case can be deduced from the power series description of the type curve:
$$ W(u) = -0.5773 - \ln(u) + u - \frac {u^2} {2 . 2!} + \frac {u^3} {3 . 3!} - \frac {u^4} {4 . 4!} + ... $$It is clear that all terms to the right of u will be smaller than u when $u<1$. Hence when u is so small that it can be neglected relative to $\ln(u)$, then all the terms to the right of $\ln(u)$ can be neglected. Therefore we have the following spproximation
$$ W(u) \approx -0.5772 -\ln(u) + O(u) $$for
$$ -\ln(u)>>u \,\,\,\rightarrow \,\,\, \ln(u)<<-u \,\,\, \rightarrow \,\,\, u<<e^{-u} \, \approx \,1 $$which is practically the case for $u<0.01$, as can also be seen in the graph for $1/u = 10^2 $. From the graph one may conclude that even for 1/u>10 or u<0.1, the logarithmic type curve is straight and therefore can be accurately computed using a logarithmic approximation of the type curve.
Below the error between the full Theis curve $W(u)$ and the approximation $Wu(u) = -0.5772 - \ln(u)$ are computed and shown. This reveals that at $u=0.01$ the error is 5.4% and at $u=0.001$ it has come down to only 0.2%.
In [56]:
U = np.logspace(-2, 0, 21)
Wa = lambda u : -0.5772 - np.log(u)
print("{:>12s} {:>12s} {:>12s} {:>12s}".format('u','W(u)','Wa(u)','1-Wa(u)/W(u)'))
print("{:>12s} {:>12s} {:>12s} {:>12s}".format(' ',' ',' ','the error'))
for u in U:
print("{:12.3g} {:12.3g} {:12.3g} {:12.1%}".format(u, W(u), Wa(u), 1-Wa(u)/W(u)))
In [57]:
U = np.logspace(-7, 1, 81)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set(xlabel='1/u', ylabel='W(u)', title='Theis type curve and its logarithmic approximation', yscale='linear', xscale='log')
ax.grid(True)
ax.plot(1/U, W(U), 'b', linewidth = 2., label='Theis type curve')
ax.plot(1/U, Wa(U), 'r', linewidth = 0.25, label='log approximation')
ax.invert_yaxis()
plt.legend(loc='best')
plt.show()
Hence, in any practical situation, the logarithmic approximation is accurate enough when $u<0.01$.
The approximatin of the Theis type curve can no be elaborated:
$$ Wa (u) \approx -0.5772 - \ln(u) = \ln(e^{-0.5772}) - \ln(u) = \ln(0.5615) - \ln(u) = \ln \frac {0.5615} {u} $$Because $u = \frac {r^2 S} {4 kD t}$ we have, with 4\times 0.5615 \approx 2.25
$$ W(u) \approx \ln \frac {2.25 kD t} {r^2 S} $$and so the drawdown approximation becomes
$$ s \approx \frac Q {4 \pi kD} \ln \frac {2.25 kD t} {r^2 S} $$The condition u<0.1 can be translated to $\frac {r^2 S} {4 kD t} < 0.1$ or
$$\frac t {r^2} > 2.5 \frac {S} {kD}$$The previous logarithmic drawdown type curve versus $1/u$ can be seen an image of the drawdown for a fixed distance and varying time. This is because $1/u$ is proportional to the real time. On the other hand, the drawdown type curve versus u may be regarded as the drawdown at a fixed time for varying distance. This follows from s versus u is
$$ W(u)\approx \ln \frac {2.25 kD t} { r^2 S} \,\,\,\, versus\,\,\,\, u = \ln \frac {r^2 S} {4 kD t} = 2 \ln \left( \frac {S} {4 kD t} r\right) $$That is, proportional r on log scale. The plot reveals this:
In [58]:
U = np.logspace(-7, 1, 81)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set(xlabel='u', ylabel='W(u)', title='Theis type curve and its logarithmic approximation', yscale='linear', xscale='log')
ax.grid(True)
ax.plot(U, W(U), 'b', linewidth = 2., label='Theis type curve')
ax.plot(U, Wa(U), 'r', linewidth = 0.25, label='log approximation')
ax.invert_yaxis()
plt.legend(loc='best')
plt.show()
This shows that the radius of influence is limited. We can now approximate this radius of influence by saying that the radius is where the appoximated Theis curve, that is the straight red line in the graph intersects the zero drawdown, i.e. $W(u) = 0$.
Hence, for the radius of influence, R, we have
$$ \ln \frac {2.25 kD t} {R^2 S} = 0 $$impying that
$$ \frac {2.25 kD t } { R^2 S } = 1 $$$$ R =\sqrt { \frac {2.25 kD t} D} $$with R the radius of influence. Computing the radius of influence is an easy way to determine how far out the drawdown affects the groundwater heads.
Below are the data given that were obtained from a pumping test carried out on the site "Oude Korendijk" south of Rotterdam in the Netherlands (See Kruseman and De Ridder, p56, 59). The piezometers are all open at 20 m below ground surface. The groundwater head is shallow, within a m from ground surface. The first18 m below ground surface consist of clay,peat and clayey fine sand. These layers form a practially impermeable confining unit. Below this, between 18 and25 m below ground surface are 7 m of sand an some gravel, that form the aquifer. Fine sandy and clayey sediments thereunder from the base of the aquifer, which is considered impermeable. Piezometers wer installed at 30, 90 and 215 m from the well, open at 20 m below ground surface. The well has its screen installed over the whole thickness of the aquifer. We consider the aquifer as confined with no leakage. But we should look with a critical eye that the drawdown curves to verify to what extent this assumption holds true.
The drawdown data for the three piezometers is given below. The first column is time after the start of the pump in minutes; the second column is the drawdown in m.
The well extracts 788 m3/d
The objective of the pumping test is to determine the properties kD and S of the aquifer.
In [59]:
# t[min], s[m]
H30 = [ [0.0, 0.0],
[0.1, 0.04],
[0.25, 0.08],
[0.50, 0.13],
[0.70, 0.18],
[1.00, 0.23],
[1.40, 0.28],
[1.90, 0.33],
[2.33, 0.36],
[2.80, 0.39],
[3.36, 0.42],
[4.00, 0.45],
[5.35, 0.50],
[6.80, 0.54],
[8.30, 0.57],
[8.70, 0.58],
[10.0, 0.60],
[13.1, 0.64]]
In [60]:
# t[min], s[m]
H90= [[0.0, 0.0],
[1.5, 0.015],
[2.0, 0.021],
[2.16, 0.23],
[2.66, 0.044],
[3.00, 0.054],
[3.50, 0.075],
[4.00, 0.090],
[4.33, 0.104],
[5.50, 0.133],
[6.0, 0.154],
[7.5, 0.178],
[9.0, 0.206],
[13.0, 0.250],
[15.0, 0.275],
[18.0, 0.305],
[25.0, 0.348],
[30.0, 0.364]]
In [61]:
# t[min], s[m]
H215=[[0.0, 0.0],
[66.0, 0.089],
[127., 0.138],
[185., 0.165],
[251., 0.186]]
Show the drawdown data on half-log time scale for the three piezometers.
What do you expect these curves to look like? Do the drawdown lines become parallel after an initial time?
Use the simplified drawdown formula to interpret the test.
Show the drawdown $s$ versus time on a double log graph for all three piezometers
Show the drawdown $s$ versus $ t/r^2 $, also for all three piezometers.
What is the difference between the last two graphs (versus $t/r^2$ instead of versus $t$) ?
Match the computed drawdown by plotting it on the same graph and adapting the transmissivity and the storage coefficient.
For your analysis, write the computed drawdown as follows:
$$ s = A\times W(u \times B) $$Then adapting A willshift the entire graph vertically, while adaptin B will shift it horizontally. This makes it easy to lay the computed curve on the data points.
By adapting A and B determine their numerical values.
Determine the transmissivity kD and the storage coefficient S from A and B.
Hint: Having determined A and B compute kD and S by matching this formula with the true Theis drawdown which is
$$ s = \frac Q {4 \pi kD} W(\frac {r^2 S} { 4 kD t}) $$That is, to make both equations equal, set $$ A = \frac Q {4 \pi kD} $$ and set $$ \frac 1 {u\, B} = \frac {4 kD} S \frac t {r^2} $$
1. How much can we pump until the head drops to halfway the screen in 3 years?
1. How much if it were 30 years?
1. How much could we pump over the same period had we drilled three of the same wells on one line at 250 m mutual distance?