9.5. Wells in a two-aquifer system

The solution for steady flow to a well in the bottom aquifer of a two-aquifer system was presented in Section 6.5. In this section, the transient equivalent is presented. The solution is used to investigate the transient head response in a two-aquifer system. When the top aquifer is unconfined, the head response in the top aquifer will be significantly slower than the head response in the bottom aquifer, also referred to as the delayed response of the water table. The solution is quite complicated as it involves solving a system of two coupled differential equations using the Laplace transform method. Although the equations are a bit longer, the Python script of the solution is relatively short.
      Consider transient flow in a two-aquifer system. The top aquifer has index 0, and the bottom aquifer has index 1. The transmissivity of each aquifer is again approximated as constant in both space and time. The storage coefficients of the top and bottom aquifers are \(S_0\) and \(S_1\), respectively. The two aquifers are separated by a leaky layer with resistance \(c\) (Figure 6.14). The storage in the leaky layer is neglected. The transient equivalent of the system of two linked differential equations ((6.63)) is

(9.52)\[ \begin{align}\begin{aligned} & \nabla^2 h_0 = \frac{S_0}{T_0}\frac{\partial h_0}{\partial t} + \frac{h_0-h_1}{cT_0}\\& \nabla^2 h_1 = \frac{S_1}{T_1}\frac{\partial h_1}{\partial t} + \frac{h_1-h_0}{cT_1} \end{aligned}\end{align} \]

A well starts pumping with a discharge \(Q\) at time \(t=0\) in aquifer 1 (the bottom aquifer). The boundary conditions are

(9.53)\[ \begin{align}\begin{aligned} \lim_{r \to 0} -2\pi r Q_{r_0} & = 0 \\\lim_{r \to 0} -2\pi r Q_{r_1} & = Q {\hskip 2em} t\ge 0 \end{aligned}\end{align} \]

(9.54)\[ h_0|_{r=\infty,t} = h_1|_{r=\infty,t} = 0 \]

and the initial condition is

(9.55)\[ h_0|_{r,t=0} = h_1|_{r,t=0} = 0 \]

A solution is obtained using Laplace transforms, as explained in Section 5.4. The solution procedure is a bit more complicated, however, as the system of two linked differential equations needs to be broken down into the solution of two unlinked differential equations. This unlinking is more difficult than for steady flow (Section 6.4) and requires the solution of an eigenvalue problem.
      The Laplace-transform of the system of differential equations (9.52) is

(9.56)\[\begin{split} \nabla^2 \overline{h}_0 = \frac{S_0 p}{T_0}\overline{h}_0 + \frac{\overline{h}_0-\overline{h}_1}{cT_0} \\ \nabla^2 \overline{h}_1 = \frac{S_1 p}{T_1}\overline{h}_1 + \frac{\overline{h}_1-\overline{h}_0}{cT_1} \end{split}\]

which may be written in matrix form as

(9.57)\[\begin{split} \nabla^2\left( \begin{array}{c} \overline{h}_0 \\ \overline{h}_1 \end{array} \right)= \left( \begin{array}{cc} \frac{S_0 c p + 1}{cT_0} & \frac{-1}{cT_0} \\ \frac{-1}{cT_1} & \frac{S_1 c p + 1}{cT_1} \end{array} \right) \left( \begin{array}{c} \overline{h}_0 \\ \overline{h}_1 \end{array} \right) \end{split}\]

The matrix on the right-hand side of the matrix equation is called the system matrix \(\mathbf{A}\) and is written as

(9.58)\[\begin{split} \mathbf{A}= \left( \begin{array}{cc} \frac{S_0 c p + 1}{cT_0} & \frac{-1}{cT_0} \\ \frac{-1}{cT_1} & \frac{S_1 c p + 1}{cT_1} \end{array} \right) = \mathbf{V} \mathbf{W} \mathbf{V}^{-1} \end{split}\]

where \(\mathbf{W}\) is a diagonal matrix with the eigenvalues of matrix \(\mathbf{A}\) on the diagonal

(9.59)\[\begin{split} \mathbf{W} = \left( \begin{array}{cc} w_0 & 0 \\ 0 & w_1 \end{array} \right) \end{split}\]

and \(\mathbf{V}\) is a matrix with the corresponding eigenvectors \(\vec{u}=\left(\begin{array}{c} u_0 \\ u_1 \end{array}\right)\) and \(\vec{v}=\left(\begin{array}{c} v_0 \\ v_1 \end{array}\right)\) as its columns

(9.60)\[\begin{split} \mathbf{V} = \left( \begin{array}{cc} u_0 & v_0 \\ u_1 & v_1 \end{array} \right) \end{split}\]

The boundary conditions in the Laplace domain are

(9.61)\[\begin{split} \begin{split} \lim_{r \to 0} -2\pi r \overline{Q}_{r_0} &= 0 \\ \lim_{r \to 0} -2\pi r \overline{Q}_{r_1} &= \frac{Q}{p} \end{split} \end{split}\]

(9.62)\[ \overline{h}_0|_{r=\infty} = \overline{h}_1|_{r=\infty} = 0 \]

The solution in the Laplace domain is

(9.63)\[\begin{split} \left( \begin{array}{c} \overline{h}_0 \\ \overline{h}_1 \end{array} \right) = -\frac{Q}{2 \pi T_1 p}\left[ b_0\text{K}_0\left(r\sqrt{w_0}\right) \left( \begin{array}{c} u_0 \\ u_1 \end{array} \right) + b_1\text{K}_0\left(r\sqrt{w_1}\right) \left( \begin{array}{c} v_0 \\ v_1 \end{array} \right)\right] \end{split}\]

where \(b_0\) and \(b_1\) are the solution of the following system of two linear equations

(9.64)\[\begin{split} \left( \begin{array}{cc} u_0 & v_0 \\ u_1 & v_1 \end{array} \right) \left( \begin{array}{c} b_0 \\ b_1 \end{array}\right) = \left( \begin{array}{c} 0 \\ 1 \end{array}\right) \end{split}\]

The Python implementation of the solution makes use of the eig function of the lingalg subpackage of numpy. The eig function returns the array \([w_0, w_1]\) of eigenvalues and the two-dimensional array \(\mathbf{V}\) with the eigenvectors.
      In the example below, the top aquifer is unconfined so that the storage coefficient \(S_0\) represents the phreatic storage of the top aquifer. The head is computed and plotted in an observation well near the well (Figure 9.10). The head responds very quickly in the bottom aquifer, but it takes approximately one day before the head in the unconfined top aquifer starts to respond. The head in the bottom aquifer starts to go down again when the head in the top aquifer starts to go down. This delayed response of the water table can be seen especially well in a semi-log plot. For early times (in this case, for the first day), the head response is similar to a Hantush well, as the head in the top aquifer varies very little at early time.

# parameters
T0 = 100 # transmissivity top aquifer, m^2/d
S0 = 0.1 # specific yield top aquifer, -
T1 = 200 # transmissivity bottom aquifer, m^2/d
S1 = 0.001 # storage coefficient bottom aquifer, -
c = 500 # resistance of leaky layer, d
Q = 800 # discharge of well, m^3/d
rw = 0.3 # radius of well, m

# solution
from scipy.special import k0
def hmaqbar(plist, r, layer):
    hbar = np.zeros(len(plist))
    for i, p in enumerate(plist):
        A = np.array([[(p * c * S0 + 1) / (c * T0), -1 / (c * T0)],
                      [-1 / (c * T1), (p * c * S1 + 1) / (c * T1)]])
        w, v = np.linalg.eig(A)
        b = np.linalg.solve(v, [0, 1])
        hbar[i] = -Q / (2 * np.pi * T1 * p) * (
            b[0] * v[layer, 0] * k0(r * np.sqrt(w[0])) + 
            b[1] * v[layer, 1] * k0(r * np.sqrt(w[1])))
    return hbar

r = 250
t = np.logspace(-3, 2, 100)
h0 = stehfest(r, t, hmaqbar, layer=0) # stehfest function from Section 5.5
h1 = stehfest(r, t, hmaqbar, layer=1)
h1_hantush = hantushvec(r, t, T1, S1, c, Q) # hantushvec from Section 9.3

# basic plot semilog
plt.semilogx(t, h0[0])
plt.semilogx(t, h1[0])
plt.semilogx(t, h1_hantush);

Figure 9.10 Drawdown in an observation well near a pumping well in the bottom aquifer of a two-aquifer system. The top aquifer is unconfined; linear time axis (left) and log-scaled time axis (right).