Investigating the character of the Theis well function

Introduction

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)))


      u            W(u)     W1(u)  
   4.0e+00   3.7794e-03   3.7793e-03
   4.0e-01   7.0238e-01   7.0238e-01
   4.0e-02   2.6813e+00   2.6812e+00
   4.0e-03   4.9482e+00   4.9482e+00
   4.0e-04   7.2472e+00   7.2471e+00
   4.0e-05   9.5495e+00   9.5493e+00
   4.0e-06   1.1852e+01   1.1852e+01
   4.0e-07   1.4155e+01   1.4154e+01
   4.0e-08   1.6457e+01   1.6456e+01
   4.0e-09   1.8760e+01   1.8759e+01
   4.0e-10   2.1062e+01   2.1061e+01

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

Theis well function as a power series

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)))


      u         W(u)       W1(u)       W2(u)  
   4.0e+00  3.7794e-03   3.7793e-03  3.7793e-03
   4.0e-01  7.0238e-01   7.0238e-01  7.0238e-01
   4.0e-02  2.6813e+00   2.6812e+00  2.6812e+00
   4.0e-03  4.9482e+00   4.9482e+00  4.9482e+00
   4.0e-04  7.2472e+00   7.2471e+00  7.2471e+00
   4.0e-05  9.5495e+00   9.5493e+00  9.5493e+00
   4.0e-06  1.1852e+01   1.1852e+01  1.1852e+01
   4.0e-07  1.4155e+01   1.4154e+01  1.4154e+01
   4.0e-08  1.6457e+01   1.6456e+01  1.6456e+01
   4.0e-09  1.8760e+01   1.8759e+01  1.8759e+01
   4.0e-10  2.1062e+01   2.1061e+01  2.1061e+01

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()


Logarithmic approximaion of the Theis type curve

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)))


           u          W(u)          Wa(u)   1-Wa(u)/W(u)
                                               the error
        0.01          4.04          4.03          0.2%
      0.0126          3.81           3.8          0.3%
      0.0158          3.58          3.57          0.4%
        0.02          3.36          3.34          0.6%
      0.0251          3.13          3.11          0.8%
      0.0316          2.91          2.88          1.1%
      0.0398          2.69          2.65          1.5%
      0.0501          2.47          2.42          2.0%
      0.0631          2.25          2.19          2.8%
      0.0794          2.03          1.96          3.8%
         0.1          1.82          1.73          5.4%
       0.126          1.62           1.5          7.5%
       0.158          1.42          1.26         10.8%
         0.2          1.22          1.03         15.5%
       0.251          1.04         0.804         22.7%
       0.316         0.867         0.574         33.8%
       0.398         0.706         0.344         51.3%
       0.501         0.558         0.114         79.7%
       0.631         0.427        -0.117        127.3%
       0.794         0.314        -0.347        210.6%
           1         0.219        -0.577        363.1%

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}$$

Radius of influence

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.

Pumping test

Introduction

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.

The data:


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]]

To work out the test:

  1. Show the drawdown data on half-log time scale for the three piezometers.

  2. What do you expect these curves to look like? Do the drawdown lines become parallel after an initial time?

  3. Use the simplified drawdown formula to interpret the test.

    1. Look where the simplified drawdown curves become zero.
    2. Determine the drawdown increase per log-cycle of time. 4 From this information determine the transmissivity and the storage coefficient.
  4. Show the drawdown $s$ versus time on a double log graph for all three piezometers

  5. Show the drawdown $s$ versus t for all three piezometers
  6. Show the drawdown $s$ versus $ t/r^2 $, also for all three piezometers.

  7. What is the difference between the last two graphs (versus $t/r^2$ instead of versus $t$) ?

  8. 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.

  1. By adapting A and B determine their numerical values.

  2. 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} $$

Exercises

  1. Consider an aquifer with constant transmissivity kD = 900 m2/d, phreatic storage coefficient S with a well that extracts Q = 2400 m3/d. Use distances between 1 and 1000 m Use times between 0.01 and 100 days
  2. Show the drawdown as a function of time for different distances where the drawdown is on linear scale and the time at logarithmic scale. Plot both the full Theis drawdown and the approximation.
  3. Show the drawdown as a function of distance for different times where the drawdown is at linear scale and the distance is at logarithmic scale. Plot both the full Theis drawdown and the approximation.
  4. Also show the radius of influence for the different times on the drawdown versus distance curve. Just show them as red dots at s=0 and r the radius of influence.
  5. A well is installed in a desert are where there are no fixed head boundaries. The well is 60 m deep and the aquifer is 80 m thick. The top of the well screen is at 40 m below ground surface and the phreatic water level at 15 m. Test pumping has revealed that the $kD = 600$ m2/d and the specific yield $Sy = 0.25$. As a first approximation ignore that the fact that the aquifer thickness gets less because of the falling water table caused by pumping.
     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?
  6. At an historic site the excavations need to be carried out in the dry. The water levels have gone up over the last 40 years due to a large reservoir that had been installed in the river by building a dam downstream. The water level needs to be lowerd by 10 m over a square are with sides of 30 m. The transmissivity is around 1200 m2/d according to test pumping and the specific yield about $Sy = 24% $.
    1. Determine the pump capacity if 4 wells are used in the corners of the exacvation which needs to be dry only two weeks after they are installed.
    2. The excavations will last for 1 year. What will be the required pump capacity after this year?