<jupyter_start><jupyter_text>Response function comparison*Developed by <NAME> and <NAME>*The purpose of this notebook is to compare several of the response functions available in Pastas.The Gamma and Hantush response function are compared with the Four Parameter function.  Gamma Step ResponseThe Gamma step response function is defined as:$$ s(t) = A  \dfrac{1}{\Gamma(n)} \int_{0}^{t} \tau^{n-1} \cdot e^{-\tau/a} d\tau $$ Hantush Step ResponseThe Hantush step response function is defined as:$$ s(t) = A  \dfrac{1}{\int_{0}^{\infty} \tau^{-1} \cdot e^{-\tau/a - b/\tau} d\tau} \int_{0}^{t} \tau^{-1} \cdot e^{-\tau/a - b/\tau} d\tau $$ FourParam Step ResponseBoth these response functions are compared with the Four Parameter response function. The Four Parameter response function is defined as:$$ s(t) = A  \dfrac{1}{\int_{0}^{\infty} \tau^{n-1} \cdot e^{-\tau/a - b/\tau} d\tau} \int_{0}^{t} \tau^{n-1} \cdot e^{-\tau/a - b/\tau} d\tau $$The Hantush and Gamma response function are special cases of the Four Parameter function. In this notebook these response functions are compared. This is done using syntheticly generated groundwater observations.First the required packages are imported and the data is loaded. The rainfall and evaporation data are imported from KNMI station De Bilt using Pastas. An extraction time series is imported using the Pandas package. For this example the time series are selected from the year 1980 until 2000.<jupyter_code>import numpy as np
from scipy.special import gammainc, gammaincinv
from scipy.integrate import quad
import pandas as pd
import pastas as ps
import matplotlib.pyplot as plt

%matplotlib inline

ps.show_versions()
rain = ps.read.read_knmi('../data/etmgeg_260.txt', variables='RH').series
evap = ps.read.read_knmi('../data/etmgeg_260.txt', variables='EV24').series
extraction = pd.read_csv('../data/nb9_extraction.csv', squeeze=True,
                         index_col=0, parse_dates=True, dayfirst=True)<jupyter_output>INFO: Inferred frequency from time series RH 260: freq=D 
INFO: Inferred frequency from time series EV24 260: freq=D<jupyter_text>Defining required functionsSeveral function are defined to generate the synthetic groundwater observations. In this example two groundwater series are generated, one using a Gamma response function and one using a Hantush response function.<jupyter_code>def gamma_tmax(A, n, a, cutoff=0.999):
    return gammaincinv(n, cutoff) * a

def gamma_step(A, n, a, cutoff=0.999):
    tmax = gamma_tmax(A, n, a, cutoff)
    t = np.arange(0, tmax, 1)
    s = A * gammainc(n, t / a)
    return s

def gamma_block(A, n, a, cutoff=0.999):
    # returns the gamma block response starting at t=0 with intervals of delt = 1
    s = gamma_step(A, n, a, cutoff)
    return np.append(s[0], s[1:] - s[:-1])

def hantush_func(t, a, b):
    return (t ** -1) * np.exp(-(a / t) - (t / b))

def hantush_step(A, a, b, tmax=1000, cutoff=0.999):
    t = np.arange(0, tmax)
    f = np.zeros(tmax)
    for i in range(1,tmax):
        f[i] = quad(hantush_func, i-1, i, args=(a, b))[0]
    F = np.cumsum(f)
    return (A / quad(hantush_func, 0, np.inf, args=(a, b))[0]) * F

def hantush_block(A, a, b, tmax=1000, cutoff=0.999):
    s = hantush_step(A, a, b, tmax=tmax, cutoff=cutoff)
    return s[1:] - s[:-1]<jupyter_output><empty_output><jupyter_text>Comparing the Gamma and the Four Parameter response functionThe first test is to compare the Gamma with the Four Parameter response function. Using the function defined above the Gamma block response function can be generated. The parameters for the block response `Atrue`, `ntrue` and `atrue` are defined together with the `dtrue` parameter. A synthetic groundwater head series is generated using the block response function and the rainfall data series.<jupyter_code>Atrue = 800
ntrue = 1.1
atrue = 200
dtrue = 20
h = gamma_block(Atrue, ntrue, atrue) * 0.001
tmax = gamma_tmax(Atrue, ntrue, atrue)
plt.plot(h)
plt.xlabel('Time (days)')
plt.ylabel('Head response (m) due to 1 mm of rain in day 1')
plt.title('Gamma block response with tmax=' + str(int(tmax)));
step = gamma_block(Atrue, ntrue, atrue)[1:]
lenstep = len(step)
h = dtrue * np.ones(len(rain) + lenstep)
for i in range(len(rain)):
    h[i:i + lenstep] += rain[i] * step
head = pd.DataFrame(index=rain.index, data=h[:len(rain)],)
head = head['1990':'1999']

plt.figure(figsize=(12,5))
plt.plot(head,'k.', label='head')
plt.legend(loc=0)
plt.ylabel('Head (m)')
plt.xlabel('Time (years)')<jupyter_output><empty_output><jupyter_text>Create Pastas modelThe synthetic head series is used as input for the Pastas model. A stress model is created with the rainfall data series and the Gamma response function. The stress model is added to the Pastas model, and the model is solved.<jupyter_code>ml = ps.Model(head)
sm = ps.StressModel(rain, ps.Gamma, name='recharge', settings='prec')
ml.add_stressmodel(sm)
ml.solve(noise=False)
ml.plots.results()<jupyter_output>INFO: Inferred frequency from time series 0: freq=D 
INFO: Inferred frequency from time series RH 260: freq=D<jupyter_text>The results of the Pastas simulation show that Pastas is able to simulate the synthetic groundwater head. The parameters calculated with Pastas are equal to the parameters used to generate the synthetic groundwater series; `Atrue`, `ntrue`, `atrue` and `dtrue`.   Create Pastas model using the Four Parameter response functionThe next step is to simulate the synthetic head series using Pastas with the Four Parameter response function.A new pastas model is created using the same head series as input. A stressmodel is created with the rainfall and the Four Parameter response function. The model is solved and the results are plotted.<jupyter_code>ml2 = ps.Model(head)
sm2 = ps.StressModel(rain, ps.FourParam, name='recharge', settings='prec')
ml2.add_stressmodel(sm2)
ml2.solve(noise=False)
ml2.plots.results()<jupyter_output>INFO: Inferred frequency from time series 0: freq=D 
INFO: Inferred frequency from time series RH 260: freq=D<jupyter_text>The results of the Pastas simulation show that the groundwater head series can be simulated using the Four Parameter resposne function. The parameters calculated using Pastas only slightly deviate from the parameters `Atrue`, `ntrue`, `atrue` and `dtrue` defined above. The parameter `recharge_b` is almost equal to 0 (meaning that the Four Parameter responce function is almost equal to the Gamma response function, as can be seen above).   Comparing the Hantush and Four Parameter response functionIn the second example of this notebook the Four Parameter response function is compared to the Hantush response function. A Hantush block response is plotted using the parameters; `Atrue_hantush`, `atrue_hantush`, `btrue_hantush`. The parameter `btrue_hantush` is calculated using `rho` and `atrue_hantush` according to Veling & Maas (2010).The Hantush block response is used together with the extraction data series to simulate the synthetic groundwater head.<jupyter_code>Atrue_hantush = -0.01 # Atrue is negative since a positive extraction results in a drop in groundwater head.
atrue_hantush = 100 # the parameter a is equal to cS in the hantush equation.
rho = 2
btrue_hantush = atrue_hantush * rho ** 2 / 4
dtrue_hantush = 20

h_hantush = hantush_block(Atrue_hantush, atrue_hantush, btrue_hantush) 
plt.plot(h_hantush)
plt.xlabel('Time (days)')
plt.ylabel('Head response (m) due to 1 m3 of extraction in day 1')
plt.title('Hantush block response with tmax=' + str(1000));
step_hantush = hantush_block(Atrue_hantush, atrue_hantush, btrue_hantush)[1:]
lenstep = len(step_hantush)
h_hantush = dtrue * np.ones(len(extraction) + lenstep)
for i in range(len(extraction)):
    h_hantush[i:i + lenstep] += extraction[i] * step_hantush
head_hantush = pd.DataFrame(index=extraction.index, data=h_hantush[:len(extraction)],)
head_hantush = head_hantush['1990':'1999']

plt.figure(figsize=(12,5))
plt.plot(head_hantush,'k.', label='head')
plt.legend(loc=0)
plt.ylabel('Head (m)')
plt.xlabel('Time (years)')<jupyter_output><empty_output><jupyter_text>Create Pastas modelA Pastas model is created using the `head_hantush` series as input. A stress model is created with the Pastas Hantush response function and the extraction as input. The stress model is added to the Pastas model and the Pastas model is solved.<jupyter_code>ml3 = ps.Model(head_hantush)
sm3 = ps.StressModel(extraction, ps.Hantush, name='extraction', 
                     settings='well', up=False)
ml3.add_stressmodel(sm3)
ml3.solve(noise=False)
ml3.plots.results()<jupyter_output>INFO: Inferred frequency from time series 0: freq=D 
INFO: Inferred frequency from time series Extraction: freq=D<jupyter_text>The results of the Pastas simulation show that the observed head can be simulated using the Hantush response function. The parameters calibrated with Pastas are very close to the true parameters.   Create Pastas model using the Four Parameter response functionA new Pastas model is created. A stress model is created using the extraction data series and the Four Parameter function, `FourParam`, as input. The stress model is added to the Pastas model an the model is solved.<jupyter_code>ml4 = ps.Model(head_hantush)
sm4 = ps.StressModel(extraction, ps.FourParam, name='extraction', 
                     settings='well', up=False)

ml4.add_stressmodel(sm4)
ml4.solve(noise=False)
ml4.plots.results()<jupyter_output>INFO: Inferred frequency from time series 0: freq=D 
INFO: Inferred frequency from time series Extraction: freq=D
{'fraction_non_alphanumeric': 0.10246719160104988, 'fraction_numerical': 0.018477690288713912, 'mean_word_length': 5.017687934301958, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '334759', 'license': 'MIT', 'n_tokens_mistral': 3206, 'n_tokens_neox': 2989, 'n_words': 1114, 'path': 'examples/notebooks/9_Response function comparison.ipynb', 'repo_name': 'RDWimmers/pastas', 'size': '9534'}