In this post, we discuss two simple models for unsteady groundwater flow in aquifers. The first is the Theis model, which is used to assess the drawdown afforded by wells in confined aquifers. The second is the Hantush-Jacob model, which applies to leaky semiconfined aquifers. The MATLAB codes for example problems 1 to 3 can be downloaded in this Google Drive folder.

**1. The Theis solution for wells in a confined aquifer**

The confined-aquifer flow equation in radial coordinates can be expressed as (see Figure 1)

where *s(r,t) *= *h _{0}* –

*h*is the drawdown,

*h*is the initial piezometric head,

_{0}*h*is the piezometric head at a radial distance

*r*from the well,

*t*is time,

*S*is the storativity (or storage coefficient) of the aquifer, and

*T*is the transmissivity of the aquifer. In 1935, C.E. Theis published a solution to equation (1) based on an analogy with radial heat transfer. The Theis model is based on several assumptions, namely:

- The aquifer is homogeneous and isotropic.
- The aquifer is uniform in thickness and the pumping never affects its exterior boundary (i.e., the aquifer can be considered infinite in areal extent).
- The aquifer is confined and does not receive any recharge.
- Well discharge is derived entirely from aquifer storage.
- The water removed from storage is discharged instantaneously when the head is lowered.
- The pumping rate is constant.
- The radius of the well is infinitely small, so that storage in the well can be ignored altogether.
- The pumping well is fully penetrating, so it receives water from the entire thickness of the aquifer and it is 100% efficient (there are no well losses).
- The initial potentiometric surface (before pumping) is horizontal.

*Figure 1. A well in a confined aquifer. *

In his derivation, Theis worked with two boundary conditions. The first BC is that the potentiometric surface formed by the well is initially horizontal,

The second BC relates to the assumption that the aquifer has an infinite horizontal extent and that there is no drawdown at very long distances from the well,

where *Q *is the pumping rate in the well. Pressing on with his solution, Theis obtained the following expression for aquifer drawdown *s*,

where *u *is given by

The integral in (4) is a modified version of the so-called exponential integral; it is a special function that appears in unsteady problems of fields such as heat transfer and neutron transport. In the mathematical literature, exponential integrals are written with the notations *Ei(u) *and *E _{1}(u)*, but groundwater scientists traditionally adopt the notation

*W(u)*, so (4) can be restated as

Equation (6) can be used to compute the drawdown afforded by a well in a confined aquifer.

**A note on the well function**

The well function can be expressed by the series expansion

which has the general form

**Example problem 1**

A well that is screened in a confined aquifer is to be pumped at a rate of 165,000 ft^{3}/day for 30 days. If the aquifer transmissivity is 5320 ft^{2}/day, and the storativity is 0.0007, what is the drawdown at distances of 50, 150, 250, 500, 1000, 3000, 5000, and 10,000 ft?

The solution to this problem is provided in code *theis.m*. As an example, we shall carry out the calculations for *r *= 500 ft. We first compute argument *u*,

The corresponding well function *W(u) *can be determined using MATLAB’s built-in function *expint*, which yields *W(u) ** *7.63. We proceed to compute the drawdown *s*,

Thus, after 30 days of pumping the drawdown at a distance of 500 ft from the well equals 18.8 ft. Similar calculations are carried out for other distances. The calculations are automated in *theis*.*m*. The drawdown-versus-distance plot is shown below.

**2. Determining aquifer transmissivity and storativity with the Cooper-Jacob method**

Importantly, the Theis solution is not just a tool for calculating well drawdowns; with the proper adjustments, it can also be used to determine the storativity and transmissivity of a confined aquifer. C.E. Jacob and H.H. Cooper observed that after a pumping well has been running for some time, *u *becomes small and the higher-power terms in the series expansion for well function *W(u) *(equation (7)) can be neglected with little loss of precision. Therefore, if

we can state that

On semilog paper, equation (10) plots as a straight line if the limiting condition is met.

In a Cooper-Jacob drawdown-time analysis, drawdown measurements are made in an observation well at various times; distance *r *is therefore constant. A plot is made between drawdown, *s *(linear scale on vertical coordinate), and time, *t* (log scale on horizontal coordinate) (Figure 2a). With reference to equation (10), the slope of the line is

or

For a change in drawdown encompassing one log cycle, we have *log(t _{2}/t_{1}) *= 1 and equation (12) reduces to

Solving for transmissivity,

In turn, where the straight line intersects the *x*-axis, the drawdown is zero and the time is *t _{0}*. Substituting these values into (10), we have

Noting that 10^{0} = 1, we can remove the logarithm and solve for storativity,

To summarize, equations (14) and (16) can be used to yield the transmissivity and storativity of a confined aquifer, respectively.

*Figure 2. Cooper-Jacob method. (a) Drawdown-time analysis. (b) Drawdown-distance analysis. (c) Drawdown-t ^{2}/r analysis. *

It is important to note that at least two variants of the Cooper-Jacob technique have been used in groundwater practice. The first variant is a *drawdown-distance *analysis, whereby drawdown measurements are made at a given time in various wells, so it is time, not distance, that is held constant. A plot is made between drawdown *s *on the vertical scale and distance *r *on a logarithmic horizontal scale (Figure 2b). From similar considerations as in drawdown-time analysis, the transmissivity is given by

and the storativity is

where *r _{0}* is the intercept at the

*x*-axis. The second – and most general – variant is to obtain measurements from several wells at various times and plot drawdown versus

*log*as shown in Figure 2c. Then, the transmissivity has the same form as equations (14) and (17), whereas the storativity is expressed as

_{10}(t/r^{2})where (*t*/*r ^{2})_{0}* is the intercept of the straight line on the

*x-*axis.

**Example problem 2**

The following data are taken from a pumping test wherein a well was pumped at a rate of 0.45 cubic feet per second. The drawdown data is for an observation well located 240 ft away from the pumped well. Plot the data on semilogarithmic paper and use the Cooper-Jacob method to estimate the transmissivity and storativity of the local aquifer.

The data are plotted using MATLAB code *cooperjacob.m*, which can be found in the Google Drive folder mentioned in the beginning of the post.

Noting that one cycle covers a drawdown interval *s *= 1.2 ft, we can use equation (14) to obtain the transmissivity *T*,

Now, the horizontal intercept of the drawdown data is *t _{0}* = 0.35 min. Noting that 1 day = 1440 min and substituting into equation (16) gives, for the storativity

*S*,

**3. The Hantush-Jacob solution for wells in a leaky, semiconfined aquifer**

Most aquifers are not isolated from sources of vertical recharge. Aquitards, either above or below the aquifer, can leak water into the aquifer if the direction of the hydraulic gradient is favorable. An aquifer overlain by a permeable confining layer is known as a *leaky aquifer*. Unsteady drawdown due to a constant-rate pumping well in this type of aquifer can be assessed with the Hantush-Jacob model. The governing equation those authors worked with was (see Figure 3)

Here, is a parameter known as the leakage factor,

where *b’ *is the thickness of the aquitard, *T *is the transmissivity of the confined aquifer, and *K’ *is the hydraulic conductivity of the aquitard. Clearly, the leakage factor is a combined property of the aquifer and aquitard. As noted by Cheng (2000), the ratio *b’*/*K’ *is an aquitard property called the *leakage coefficient *or *leakance*. Its inverse, *b’*/*K’*, is known as the *hydraulic resistance*. The leakage coefficient represents the aquitard flow across a unit horizontal area into the adjacent aquifer under a unit head difference between the top and bottom of the aquitard.

Hantush and Jacob assumed that all the water flowing to the pumping well came from either elastic storage in the confined aquifer or leakage across the confining layer. No water is derived from elastic storage in the confining layer. To summarize, the assumptions used by those authors include:

- The aquifer is bounded on the top by an aquitard.
- The aquitard is overlain by an unconfined aquifer, known as the
*source bed*. - The water table in the source bed is initially horizontal.
- The water table in the source bed does not fall during pumping of the aquifer.
- Groundwater flow in the aquitard is vertical.
- The aquifer is compressible, and water drains instantaneously with a decline in head.
- The aquitard is incompressible, so that no water is released from storage in the aquitard when the aquifer is pumped.

The solution to (20) obtained by Hantush and Jacob is

Here, *W _{HJ}* is the Hantush-Jacob leaky aquifer well function, which is expressed as

where *u *is the same parameter used in the Theis solution.

*Figure 3. Well in a leaky aquifer. *

**Example problem 3**

If the aquifer described in example problem 1 is not fully confined, but is overlain by a 8.0-ft thick leaky, confining layer with a vertical hydraulic conductivity of 0.034 ft/day, what would be the drawdown values after 30 days of pumping at 165,000 ft^{3}/day at the indicated distances?

As an example, we shall carry out the calculations for *r *= 500 ft. We first compute argument *u* (whose value is identical to the one obtained in problem 1),

Next, we compute the leakage factor ,

Then, we proceed to compute the leaky aquifer well function *W _{HJ}(u,*r/

*)*. There is no analytical solution to the integral that defines

*W*

_{HJ}(u, r/*)*, but we can easily approximate the integral numerically; in code

*hantushjacob.m,*we use MATLAB’s Symbolic Math Toolbox to declare the integral symbolically and then apply the function

*vpa()*to evaluate the integral numerically. For the example at hand,

Lastly, the drawdown *s *is determined as

Thus, after 30 days of pumping the drawdown at a distance of 500 ft from the well equals 5.03 ft. Similar calculations apply to other values of radial distance. The calculations are automated in *hantushjacob*.*m*. The drawdown-versus-distance plot is shown below.

**References**

- CHENG, A.H.-D. (2000).
*Multilayered Aquifer Systems*. Marcel Dekker. - FETTER, C.W., Jr. (2014).
*Applied Hydrogeology.*4th edition. Pearson. - GUPTA, R.S. (2017).
*Hydrology and Hydraulic Systems.*4th edition. Waveland Press.