Check Response Functions#

This notebook has two topics: 1. It checks the step response functions by numerically integrating the impulse response functions. 2. It compares the Polder, Hantush, and HantushWell response functions to their classical equivalents and gives relationships between the parameters used in Pastas and the aquifer parameters used in the classical functions.

[1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import pastas as ps
[21]:
plt.rcParams["figure.figsize"] = [3.6, 2.4]
plt.rcParams["font.size"] = 8

Comparison of step function with numerical integration of the impulse response function#

Gamma#

[22]:
ps.Gamma.impulse
[22]:
$$ \displaystyle \begin{array}{l} \left( A\space,\space n\space,\space a\right) = p \\ \mathrm{theta}(t, p) = \frac{A t^{n - {1}} \exp{\left({\frac{-t}{a}}\right)}}{a^{n} \mathrm{{\Gamma}}\left(n\right)} \end{array} $$
[23]:
A = 5
n = 1.5
a = 50
p = [A, n, a]

gamma = ps.Gamma()
tmax = gamma.get_tmax(p)
t = np.arange(0, tmax)

step = gamma.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(gamma.impulse, 0, t[i], args=(p))[0]
[24]:
plt.plot(t[1:], step, label="analytic")
plt.plot(t, stepnum, "--", label="numerical")
plt.xlabel("time (d)")
plt.ylabel("step (m)")
plt.grid()
_ = plt.legend()  # try to show figure in readthedocs
../_images/benchmarks_check_response_functions_7_0.png

Exponential#

[25]:
ps.Exponential.impulse
[25]:
$$ \displaystyle \begin{array}{l} \left( A\space,\space a\right) = p \\ \mathrm{theta}(t, p) = \frac{A}{a} \exp{\left({\frac{-t}{a}}\right)} \end{array} $$
[26]:
A = 5
a = 50
p = [A, a]

exponential = ps.Exponential()
tmax = exponential.get_tmax(p)
t = np.arange(0, tmax)

step = exponential.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(exponential.impulse, 0, t[i], args=(p))[0]
[27]:
plt.plot(t[1:], step, label="analytic")
plt.plot(t, stepnum, "--", label="numerical")
plt.xlabel("time (d)")
plt.ylabel("step (m)")
plt.grid()
_ = plt.legend()  # try to show figure in readthedocs
../_images/benchmarks_check_response_functions_11_0.png

Hantush#

[28]:
ps.Hantush.impulse
[28]:
$$ \displaystyle \begin{array}{l} \left( A\space,\space a\space,\space b\right) = p \\ \mathrm{theta}(t, p) = \frac{A}{{2} t \mathrm{K_0}\left({2} \sqrt{b}\right)} \exp{\left({\frac{-t}{a} - \frac{a b}{t}}\right)} \end{array} $$
[29]:
A = 5
a = 50
b = 2
p = [A, a, b]

hantush = ps.Hantush()
tmax = hantush.get_tmax(p)
t = np.arange(0, tmax)

step = hantush.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(hantush.impulse, 0, t[i], args=(p))[0]
[30]:
plt.plot(t[1:], step, label="analytic")
plt.plot(t, stepnum, "--", label="numerical")
plt.xlabel("time (d)")
plt.ylabel("step (m)")
plt.grid()
_ = plt.legend()  # try to show figure in readthedocs
../_images/benchmarks_check_response_functions_15_0.png

Polder#

[31]:
ps.Polder.impulse
[31]:
$$ \displaystyle \begin{array}{l} \left( A\space,\space a\space,\space b\right) = p \\ \mathrm{theta}(t, p) = A \sqrt{\frac{a b}{{\pi}}} t^{-{1.5}} \exp{\left({\frac{-t}{a} - \frac{a b}{t}}\right)} \end{array} $$
[32]:
A = 5
a = 100
b = 0.25
p = [A, a, b]

polder = ps.Polder()
tmax = polder.get_tmax(p)
t = np.arange(0, tmax)

step = polder.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(polder.impulse, 0, t[i], args=(p))[0]
[33]:
plt.plot(t[1:], step, label="analytic")
plt.plot(t, stepnum, "--", label="numerical")
plt.xlabel("time (d)")
plt.ylabel("step (m)")
plt.grid()
_ = plt.legend()  # try to show figure in readthedocs
../_images/benchmarks_check_response_functions_19_0.png

Four-parameter function#

[34]:
ps.FourParam.impulse
[34]:
$$ \displaystyle \begin{array}{l} \left( _\space,\space n\space,\space a\space,\space b\right) = p \\ \mathrm{theta}(t, p) = t^{n - {1}} \exp{\left({\frac{-t}{a} - \frac{a b}{t}}\right)} \end{array} $$
[35]:
A = 1  # impulse response implemented for A=1 only
n = 1.5
a = 50
b = 10
p = [A, n, a, b]

fourparam = ps.FourParam(quad=False)  # use simple integration
tmax = fourparam.get_tmax(p)
t = np.arange(0, tmax)

step = fourparam.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(fourparam.impulse, 0, t[i], args=(p))[0]
stepnum = (
    stepnum / quad(fourparam.impulse, 0, np.inf, args=p)[0]
)  # four param is scaled at the end
[36]:
plt.plot(t[1:], step, label="analytic")
plt.plot(t, stepnum, "--", label="numerical")
plt.xlabel("time (d)")
plt.ylabel("step (m)")
plt.grid()
_ = plt.legend()  # try to show figure in readthedocs
../_images/benchmarks_check_response_functions_23_0.png

Double exponential function#

[37]:
ps.DoubleExponential.impulse
[37]:
$$ \displaystyle \begin{array}{l} \left( A\space,\space {\alpha}\space,\space a_1\space,\space a_2\right) = p \\ \mathrm{theta}(t, p) = A \left( \frac{{1} - {\alpha}}{a_1} \exp{\left({\frac{-t}{a_1}}\right)} + \frac{{\alpha}}{a_2} \exp{\left({\frac{-t}{a_2}}\right)} \right) \end{array} $$
[38]:
A = 5  # impulse response implemented for A=1 only
a = 10
b = 50
f = 0.4
p = [A, f, a, b]

doubexp = ps.DoubleExponential()
tmax = doubexp.get_tmax(p)
t = np.arange(0, tmax)

step = doubexp.step(p)
stepnum = np.zeros(len(t))
for i in range(1, len(t)):
    stepnum[i] = quad(doubexp.impulse, 0, t[i], args=(p))[0]
[39]:
plt.plot(t[1:], step, label="analytic")
plt.plot(t, stepnum, "--", label="numerical")
plt.xlabel("time (d)")
plt.ylabel("step (m)")
plt.grid()
_ = plt.legend()  # try to show figure in readthedocs
../_images/benchmarks_check_response_functions_27_0.png

Kraijenhoff#

Kraijenhoff van de Leur#

Impulse Response#

from A study of non-steady groundwater flow with special reference to a reservoir-coefficient (1958) formula 2

$ \theta`(t) = :nbsphinx-math:frac{4N}{pi S}` \sum_{n=1,3,5…}^:nbsphinx-math:infty `:nbsphinx-math:left`( \frac{1}{n} \exp{\left( {-n^2\frac{\pi^2T}{SL^2} t} \right)} \sin `:nbsphinx-math:left`(\frac{n\pi x}{L}\right) \right) $

Step Response#

The step response is obtained by taking the integral of the impulse response function

$ \Theta`(t) = :nbsphinx-math:frac{4 N}{pi S}` \sum_{n=1,3,5…}^:nbsphinx-math:infty `:nbsphinx-math:frac{1}{n^3}` \left`(:nbsphinx-math:frac{SL^2}{pi^2 T}` - \frac{SL^2}{\pi^2 T} \exp\left`(-n^2:nbsphinx-math:frac{pi^2T}{SL^2}`t:nbsphinx-math:right):nbsphinx-math:right) \sin `:nbsphinx-math:left`(\frac{n\pi x}{L}\right) $

$ \Theta`(t) = :nbsphinx-math:frac{4 N L^2}{pi^3 T}` \sum_{n=1,3,5…}^:nbsphinx-math:infty `:nbsphinx-math:frac{1}{n^3}` \left`(1 - :nbsphinx-math:exp`:nbsphinx-math:left`(-n^2:nbsphinx-math:frac{pi^2T}{SL^2}`t:nbsphinx-math:right):nbsphinx-math:right) \sin `:nbsphinx-math:left`(\frac{n\pi x}{L}\right)$

And \(\sum_{n=1,3,5...}^\infty n = \sum_{n=0}^\infty (2n+1)\) gives:

$ \Theta`(t) = :nbsphinx-math:frac{4 N L^2}{pi^3 T}` \sum_{n=0}^:nbsphinx-math:infty `:nbsphinx-math:frac{1}{(2n+1)^3}` \left`(1 - :nbsphinx-math:exp`:nbsphinx-math:left`(-(2n+1)^2:nbsphinx-math:frac{pi^2T}{SL^2}`t):nbsphinx-math:right):nbsphinx-math:right) \sin `:nbsphinx-math:left`(\frac{(2n+1)\pi x}{L}\right)$

Kraijenhoff van de Leur takes \(\frac{x}{L}=\frac{1}{2}\) as the middle of the domain.

Bruggeman#

from Analytical Solutions of Geohydrological Problems (1999) formula 133.15

Step Response#

$ \Theta`(t) = :nbsphinx-math:frac{-N}{2T}`:nbsphinx-math:left`(x^2 - :nbsphinx-math:frac{1}{4}`L^2:nbsphinx-math:right) - \frac{4NL^2}{\pi^3T} \sum_{n=0}^:nbsphinx-math:infty \frac{(-1)^n}{(2n + 1)^3} \cos\left`(:nbsphinx-math:frac{(2n+1)pi x}{L}`:nbsphinx-math:right) \exp\left`(-:nbsphinx-math:frac{(2n+1)^2pi^2 T}{SL^2}`t:nbsphinx-math:right) $

$ \Theta`(t) = :nbsphinx-math:frac{-NL^2}{2T}`:nbsphinx-math:left`(:nbsphinx-math:left`(\frac{x}{L}\right)^2 - \frac{1}{4}\right) - \frac{4NL^2}{\pi^3T} \sum_{n=0}^:nbsphinx-math:infty \frac{(-1)^n}{(2n + 1)^3} \exp\left`(-:nbsphinx-math:frac{(2n+1)^2pi^2 T}{SL^2}`t:nbsphinx-math:right) \cos\left`(:nbsphinx-math:frac{(2n+1)pi x}{L}`:nbsphinx-math:right) $

$ \Theta`(t) = :nbsphinx-math:frac{-NL^2}{2T}`:nbsphinx-math:left`(:nbsphinx-math:left`(\frac{x}{L}\right)^2 - \tfrac{1}{4}\right) \left`(1 - :nbsphinx-math:frac{8}{pi^3 left(frac{1}{4} - left(frac{x}{L}right)^2right)}` \sum_{n=0}^:nbsphinx-math:infty \frac{(-1)^n}{(2n + 1)^3} \exp\left`(-:nbsphinx-math:frac{(2n+1)^2pi^2 T}{SL^2}`t:nbsphinx-math:right) \cos\left`(:nbsphinx-math:frac{(2n+1)pi x}{L}`:nbsphinx-math:right) \right) $

Note that \(x=0\) is the middle of the domain for Bruggeman.

Pastas Implementation#

In Pastas the Bruggeman response function is computed and the parameters are transformed to:

Scale parameter (such that the gain is always \(A\)):

\(A = \frac{-NL^2}{2T}\left(\left(\frac{x}{L}\right)^2 - \tfrac{1}{4}\right)\)

Reservoir coefficient (also known as \(j\) in Kraijenhoff):

\(a = \frac{SL^2}{\pi^2 T}\)

Location in the domain:

\(b = \frac{x}{L}\)

Such that the step response becomes:

$ \Theta`(t) = A:nbsphinx-math:left`(1 - \frac{8}{\pi^3(\frac{1}{4} - b^2)} \sum_{n=0}^:nbsphinx-math:infty `:nbsphinx-math:frac{(-1)^n}{(2n+1)^3}` \cos\left`((2n+1):nbsphinx-math:pi b:nbsphinx-math:right`):nbsphinx-math:exp\left`(-:nbsphinx-math:frac{(2n+1)^2t}{a}`:nbsphinx-math:right) \right)$

Taking the derivative gives the impulse response:

[40]:
ps.Kraijenhoff.impulse
[40]:
$$ \displaystyle \begin{array}{l} \left( A\space,\space a\space,\space b\right) = p \\ nterms = {10} \\ \mathrm{theta}(t, p) = \frac{A {8}}{{\pi}^{{3}} \left( \frac{{1}}{{4}} - b^{{2}} \right)} \sum_{n = {0}}^{{nterms - 1}} \left({\frac{\left( -{1} \right)^{n}}{a \left( {2} n + {1} \right)} \cos{\left({\left( {2} n + {1} \right) {\pi} b}\right)} \exp{\left({\frac{-\left( \left( {2} n + {1} \right)^{{2}} t \right)}{a}}\right)}}\right) \end{array} $$
[41]:
A = 5
a = 10
b = 0.25
p = [A, a, b]

khoff = ps.Kraijenhoff()
tmax = khoff.get_tmax(p)
t = np.arange(0, tmax)

step = khoff.step(p)
stepnum_brug = np.zeros(len(t))

for i in range(1, len(t)):
    stepnum_brug[i] = quad(khoff.impulse, 0, t[i], args=(p))[0]
[43]:
plt.plot(t[1:], step, label="analytic")
plt.plot(t, stepnum_brug, "--", label="numerical Bruggeman")
plt.xlabel("time (d)")
plt.ylabel("step (m)")
plt.grid()
_ = plt.legend()  # try to show figure in readthedocs
../_images/benchmarks_check_response_functions_31_0.png

Comparison to classical analytical response functions#

Polder step function compared to classic polder function#

The classic polder function is (Eq. 123.32 in Bruggeman, 1999)

\[h(t) = \Delta h \text{P}\left(\frac{x}{2\lambda}, \sqrt{\frac{t}{cS}}\right)\]

where P is the polder function.

[44]:
from scipy.special import erfc


def polder_classic(t, x, T, S, c):
    X = x / (2 * np.sqrt(T * c))
    Y = np.sqrt(t / (c * S))
    rv = 0.5 * np.exp(2 * X) * erfc(X / Y + Y) + 0.5 * np.exp(-2 * X) * erfc(X / Y - Y)
    return rv


delh = 2
T = 20
c = 5000
S = 0.01
x = 400
x / np.sqrt(c * T)
t = np.arange(1, 121)
h_polder_classic = np.zeros(len(t))
for i in range(len(t)):
    h_polder_classic[i] = delh * polder_classic(t[i], x=x, T=T, S=S, c=c)
#
A = delh
a = c * S
b = x**2 / (4 * T * c)
hpd = polder.step([A, a, b], dt=1, cutoff=0.95)
#
plt.plot(t, h_polder_classic, label="Polder classic")
plt.plot(np.arange(1, len(hpd) + 1), hpd, label="Polder Pastas", linestyle="--")
plt.legend()
[44]:
<matplotlib.legend.Legend at 0x7fa35ccbbf70>
../_images/benchmarks_check_response_functions_34_1.png

Hantush step function compared to classic Hantush function#

The classic Hantush function is

\[h(r, t) = \frac{-Q}{4\pi T}\int_u ^\infty \exp\left(-y - \frac{r^2}{4 \lambda^2 y} \right) \frac{\text{d}y}{y}\]

where

\[u=\frac{r^2 S}{4 T t}\]

For large time, the classic Hantush function goes to

\[\lim_{t\to\infty} h(r, t) = \frac{-Q}{2\pi T}\text{K}_0(r/\lambda)\]

where \(\lambda^2=cT\). The classic Hantush function is a step function. The impulse response function \(\theta\) is obtained by differentiation

\[\theta(t) = \frac{\partial h}{\partial t} = \frac{-Q}{4\pi Tt} \exp\left(\frac{-t}{cS} -\frac{r^2cS}{4cTt}\right)\]

The Hantush impulse response function used in Pastas is defined as

\[\theta(t) = A\frac{1}{2\text{K}_0(2\sqrt{b})} \frac{1}{t} e^{-t/a - ab/t}\]

Hence, the parameters in Pastas are related to the classic Hantush function as

\[A = -\frac{\text{K}_0(2\sqrt{b})}{2\pi T}\]
\[a = cS\]
\[b = \frac{r^2}{4\lambda^2}\]
[45]:
from scipy.integrate import quad
from scipy.special import k0


def integrand_hantush(y, r, lab):
    return np.exp(-y - r**2 / (4 * lab**2 * y)) / y


def hantush_classic(t=1, r=1, Q=1, T=100, S=1e-4, c=1000):
    lab = np.sqrt(T * c)
    u = r**2 * S / (4 * T * t)
    F = quad(integrand_hantush, u, np.inf, args=(r, lab))[0]
    return -Q / (4 * np.pi * T) * F


c = 1000  # d
S = 0.01  # -
T = 100  # m^2/d
r = 500  # m
Q = 20  # m^3/d
lab = np.sqrt(T * c)
#
A = k0(r / lab) / (2 * np.pi * T)
a = c * S
b = r**2 / (4 * T * c)

#
ht = hantush.step([A, a, b], dt=1, cutoff=0.99) * -Q  # multiply by Q

# calculate classic Hantush
t = np.arange(1, len(ht) + 1)
h_hantush_classic = np.zeros(len(t))
for i in range(len(t)):
    h_hantush_classic[i] = hantush_classic(t[i], r=r, Q=Q, T=T, S=S, c=c)

# plot comparison
plt.plot(t, h_hantush_classic, label="Hantush classic")
plt.plot(np.arange(1, len(ht) + 1), ht, "--", label="Hantush Pastas")
plt.axhline(-Q * A, ls="dashed", c="k", label="Q*A")
plt.grid()
plt.legend();
../_images/benchmarks_check_response_functions_36_0.png

HantushWellModel step function in Pastas compared to classic Hantush function and Hantush function in Pastas#

The impulse response of the classic Hantush function is

\[\theta(t) = \frac{\partial h}{\partial t} = \frac{-Q}{4\pi Tt} \exp\left(\frac{-t}{cS} -\frac{r^2cS}{4cTt}\right)\]

The impulse response used for the HantushWellModel is

The HantushWellModel impulse response function used in Pastas is defined as

\[\theta(r, t) = \frac{A^\prime}{2t} e^{-t/a - a r^2 \exp(b^\prime) /t}\]

Where $ A^:nbsphinx-math:prime `$, :math:`a, and \(b^\prime\) are fitting parameters. The parameters are named $A^:nbsphinx-math:prime `$ and :math:`b^prime to distinguish them from \(A\) and \(b\) in the original Hantush formulation. The fitting parameters are related to the geohydrological parameters as

\[A^\prime = -\frac{1}{2\pi T}\]
\[a = cS\]
\[b^\prime = \ln \left( \frac{1}{4\lambda^2} \right)\]

and the gain of the HantushWellModel is

\[\text{gain} = A^\prime\text{K}_0 \left( 2 r \exp \left( \frac{b^\prime}{2} \right) \right)\]
[46]:
# compare pastas.Hantush to pastas.HantushWellModel

hantush = ps.Hantush()
hantush.update_rfunc_settings(up=False)
hantushwm = ps.HantushWellModel()
hantushwm.update_rfunc_settings(up=False)
hantushwm.set_distances(1.0)

bprime = np.log(1 / (4 * T * c))
Aprime = A / k0(2 * r * np.exp(bprime / 2))

h = hantush.step([A, a, b], dt=1, cutoff=0.99)
hwm = hantushwm.step([Aprime, a, bprime, r], dt=1, cutoff=0.99)

plt.plot(hwm, "-", label="HantushWellModel")
plt.plot(h, ".", label="Hantush")
plt.axhline(A, ls="dashed", color="k", label="Gain")
plt.legend()
plt.grid();
../_images/benchmarks_check_response_functions_39_0.png