The aquifer is has semi-infinite extent, i.e. $0 \le x \le \infty$. It is considered of constant transmissivity $kD$ and constnat storage coefficient $S$.
The partial differential equation that is governing the flow reads
$$ kD \frac {\partial^2 s} {\partial x} = S \frac {\partial s} {\partial t} $$The analytical transient solution for change of head caused by a sudden, forced, change of hat at $x=0$ and $t=0$ is equal to $A [m]$, is
$$ s(x, t) = A \,\mbox{erfc}(u), \,\,\,\, u=\sqrt{\frac {x^2 S} {4 kD t}} $$Where $\mbox{erfc} () $ is the so-called complementary error function:
$$ \mbox{erfc} (z) = \frac 2 {\sqrt {\pi} } \intop _z ^\infty e ^{-y^2}dy $$And so its derivative is
$$ \frac {d \mbox{erfc}(z)} {d z} = - \frac 2 {\sqrt {\pi}} e ^{-z^2} $$Therefore, the discharge equals
$$ Q = -kD \frac {\partial s} {\partial x} = A \sqrt{\frac {kDS} {\pi t}} \exp \left( -\frac {x^2 S} {4 kD t} \right) $$and for $ x = 0 $
$$ Q_0 = A \sqrt{\frac {kD S} {\pi t}}$$
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc
In [3]:
def newfig(title='?', xlabel='?', ylabel='?', xlim=None, ylim=None,
xscale='linear', yscale='linear', size_inches=(14, 8)):
'''Setup a new axis for plotting'''
fig, ax = plt.subplots()
fig.set_size_inches(size_inches)
ax.set_title(title)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_xscale(xscale)
ax.set_yscale(yscale)
if xlim is not None: ax.set_xlim(xlim)
if ylim is not None: ax.set_ylim(ylim)
ax.grid(True)
return ax
In [4]:
# aquifer properties
kD = 900 # m2/d
S = 0.1 # [-]
A = 1.5 # m amount of sudden rise of head at x=0 and t=0
x = np.linspace(0, 1000, 101) # coordinates at which to compute the head change and discharge
times = [0.1, 1, 2, 4, 8, 16, 32, 64] # times at which to show the had and discharges as a function of x
ax = newfig("Effect of sudden head change at $x$=0 and $t$=0", 'x [m]', 'head change s [m]')
for t in times:
u = np.sqrt((x**2 * S)/(4 * kD * t))
s = A * erfc(u)
ax.plot(x, s, label=f't = {t:.0f} d')
ax.legend()
# plot discharges
ax = newfig("Discharge due to sudden head change at $x$=0 and $t$=0", 'x [m]', 'Discharge Q [m2/d]')
for t in times:
u = np.sqrt((x**2 * S)/(4 * kD * t))
Q = A * np.sqrt((kD * S)/(np.pi * t)) * np.exp(-u**2)
ax.plot(x, Q, label='t = {:.0f} d'.format(t))
ax.legend()
Out[4]:
In [ ]: