Symmetric solution of a decaying head in strip of land

IHE, Delft, 2019-01-02

@T.N.Olsthoorn, 2019-01-02

A solution, which shows the deline of the head in a strip due to bleeding to the fixed heads at both ends.

$$ s(x, t) = A \frac 4 \pi \sum _{j=1} ^\infty \left\{ \frac {(-1)^{j-1}} {2 j - 1} \cos \left[(2 j -1) \frac \pi 2 \frac x b \right] \exp \left[ -(2 j - 1)^2 \left( \frac \pi 2 \right)^2 \frac {kD} {b^2 S} t \right] \right\} $$

To make sure that both reflex decay of an initial head of $a$ above the fixed heads at $x \pm L/2 = b$, we have to subtract the sudden head solutions from the initial head $a$



In [3]:

    
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

Convenience function for setting up graphs



In [4]:

    
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 [5]:

    
# Aquifer
kD = 600 # m2/d
S = 0.1 # [-]

# Size of strip of land
L = 150 # m (strip width)
b =L/2 # [m] half width of strip

# Chosing coordinates inside the strip of land
x = np.linspace(-L/2, L/2, 201) # points, taking left at zero.

a = 1.0 # m, sudden head change at x = -L/2
times = np.linspace(0, 0.5, 11) # 
times[0] = 1e-6 # prevent division by zero



In [8]:

    
ax = newfig('Symmetric solution for head decay in strip', 'x [m]', 'head [h]', xlim=(-b, b))

for t in times:
    h = np.zeros_like(x)
    for j in range(1,20):
        h += a * 4  / np.pi * ((-1)**(j-1) / (2 * j - 1) *
            np.cos((2 * j - 1) * np.pi / 2 * x / b) *
            np.exp(- (2 * j - 1)**2 * (np.pi / 2)**2 * kD /(b**2 * S) * t))
    ax.plot(x, h, label='t={:.1f}'.format(t))
ax.legend()









    Out[8]:





<matplotlib.legend.Legend at 0x109322e10>

As can be seen, we did not compute enough terms in the series to remove the waves in the first graph. It can be solved by taking more terms or by choosing a somewhat later first time. Just change the value in the loop (20) to for instance 200 to see the effect.

The individual terms of the analytic solution

Let's take a time and show how the individual terms add up to yield the final solution.



In [19]:

    
t = 0.02 # d

ax = newfig(f'Individual terms at t = {t:.1f} d', 'x [m]', 'head [h]', xlim=(-b, b))

t = 0.001

H = np.zeros_like(x)
for j in range(1,10):
    h = a * 4  / np.pi * ((-1)**(j-1) / (2 * j - 1) *
        np.cos((2 * j - 1) * np.pi / 2 * x / b) *
        np.exp(- (2 * j - 1)**2 * (np.pi / 2)**2 * kD /(b**2 * S) * t))
    ax.plot(x, h, label=f'term[{j:}]')
    H += h
ax.plot(x, H, lw=3, label=f'sum, t={t:.2f} d')
ax.legend()









    Out[19]:





<matplotlib.legend.Legend at 0x81d2c79d0>

This example shows how the final result is being built up from the consines, that are at ever higher frequency and ever larger damping. For larger values of $t$, only the first term will still matter. It can also be seen that the sum of the 9 terms is still wavy, so that more terms are needed to get the true result at time $t$.

At the left and right edgte, all cosines add up to yield the steep head near the boundaries of the strip.



In [ ]: