Sudden change of head at $x=0$, 1D flow in a semi-infinite aquifer

IHE, Delft, 2017

@T.N.Olsthoorn

Setting and analytical solution

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

Importing modules


In [2]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

Convenience function for setting up a graph


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

Implementation


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]:
<matplotlib.legend.Legend at 0x818cb1050>

In [ ]: