# Changing response functions#

R.A. Collenteur, University of Graz, 2021

In this notebook the new ChangeModel is tested, based on the work by Obergfjell et al. (2019). The main idea is to apply different response functions for two different periods. As an example we look at the the groundwater levels measured near the river the Mur in Austria, where a dam was recently built.

[1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pastas as ps

ps.set_log_level("ERROR")
ps.show_versions()

Python version: 3.10.8
NumPy version: 1.23.5
Pandas version: 2.0.1
SciPy version: 1.10.1
Matplotlib version: 3.7.1
Numba version: 0.57.0
LMfit version: 1.2.1
Latexify version: Not Installed
Pastas version: 1.0.1


[2]:

prec = pd.read_csv("data_step/prec.csv", index_col=0, parse_dates=True).squeeze()
river -= river.min()

axes = ps.plots.series(
stresses=[prec, evap, river],
tmin="2000",
)


# 2. The weighting factor#

The stress is convoluted two times with different response functions. Then, a weighting function is used to add the two contributions together and compute the final contribution.

[3]:

npoints = 100

tchange = 50 / npoints
t = np.linspace(0, 1, npoints)
color = plt.cm.viridis(np.linspace(0, 1, 10))

for beta, c in zip(np.linspace(-1, 1, 10), color):
beta1 = beta * npoints
omega = 1 / (np.exp(beta1 * (t - tchange)) + 1)
plt.plot(omega, color=c, label="$beta$={:.2f}".format(beta))

plt.ylabel("$\omega$ [-]")
plt.xlabel("Time [t]")
plt.legend()

[3]:

<matplotlib.legend.Legend at 0x7fbfbe7a26e0>


# 3. Make a model#

We now make two models:

• one model where we assume the response of the heads to the river level remains the same

• and one model where the response to the river levels changes.

[4]:

# Normal Model

sm = ps.StressModel(river, ps.Exponential(), name="test")
step = ps.StepModel("2012-01-01", rfunc=ps.One(), name="step")

ml.solve(report=False, tmin="2004", tmax="2017-12-31", noise=True)
ml.plots.results(figsize=(10, 6))

# ChangeModel

cm = ps.ChangeModel(
river,
ps.Exponential(),
ps.Exponential(),
name="test",
tchange="2012-01-01",
)

ml2.solve(report=False, tmin="2004", tmax="2017-12-31", noise=True)
ml2.plots.results(figsize=(10, 6));

/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.10/site-packages/pastas/stressmodels.py:1687: RuntimeWarning: overflow encountered in exp
omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.10/site-packages/pastas/stressmodels.py:1687: RuntimeWarning: overflow encountered in exp
omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.10/site-packages/pastas/stressmodels.py:1687: RuntimeWarning: overflow encountered in exp
omega = 1 / (np.exp(beta * (t - sigma)) + 1)


The second model shows a better fit, but also the step trend changed.

[5]:

print("RMSE for the first model:", ml.stats.rmse().round(2))
print("RMSE for the second model:", ml2.stats.rmse().round(2))

RMSE for the first model: 0.17
RMSE for the second model: 0.13

/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.10/site-packages/pastas/stressmodels.py:1687: RuntimeWarning: overflow encountered in exp
omega = 1 / (np.exp(beta * (t - sigma)) + 1)

[6]:

ml2.parameters

[6]:

initial name optimal pmin pmax vary stderr
test_1_A 1.00000 test 0.881948 1.000000e-05 100.0 True 9.400205e-03
test_1_a 10.00000 test 1.272077 1.000000e-02 1000.0 True 4.395089e-02
test_2_A 1.00000 test 0.428230 1.000000e-05 100.0 True 1.161800e-02
test_2_a 10.00000 test 1.930610 1.000000e-02 1000.0 True 1.026581e-01
test_beta 0.00000 test 781566.435302 -inf inf True 1.394655e-16
test_tchange 734503.00000 test 734503.000000 6.124110e+05 825914.0 False NaN
step_d 1.00000 step 1.130024 0.000000e+00 NaN True 1.750598e-02
step_tstart 734503.00000 step 734503.000000 6.124110e+05 825914.0 False NaN
constant_d 330.18797 constant 329.256815 NaN NaN True 1.046840e-02
noise_alpha 7.00000 noise 24.476228 1.000000e-05 5000.0 True 1.370048e+00

# 4. Compare the response functions#

We can also look at the response to the river before and after,

[7]:

cm_rf1 = cm.rfunc1.step(
p=ml2.parameters.loc[["test_1_A", "test_1_a"], "optimal"].values
)
cm_rf2 = cm.rfunc2.step(
p=ml2.parameters.loc[["test_2_A", "test_2_a"], "optimal"].values
)

plt.plot(np.arange(1, len(cm_rf1) + 1), cm_rf1)
plt.plot(np.arange(1, len(cm_rf2) + 1), cm_rf2)

plt.legend(["Before", "After"]);


# 5. Another way#

We can also add the stress twice, saving one parameter that needs to be estimated.

[8]:

ml3 = ps.Model(head, name="linear")

river1 = river.copy()
river1.loc["2012":] = 0

river2 = river.copy()
river2.loc[:"2011"] = 0

r1 = ps.StressModel(river1, rfunc=ps.Exponential(), name="river")
r2 = ps.StressModel(river2, rfunc=ps.Exponential(), name="river2")
step = ps.StepModel("2012-01-01", rfunc=ps.One(), name="step")

ml3.solve(report=False, tmin="2004", tmax="2017-12-31", noise=True)
ml3.plots.results(figsize=(10, 6));


# How do the results compare?#

[9]:

# change model
cm_rf1 = cm.rfunc1.step(
p=ml2.parameters.loc[["test_1_A", "test_1_a"], "optimal"].values
)
cm_rf2 = cm.rfunc2.step(
p=ml2.parameters.loc[["test_2_A", "test_2_a"], "optimal"].values
)
plt.plot(np.arange(1, len(cm_rf1) + 1), cm_rf1)
plt.plot(np.arange(1, len(cm_rf2) + 1), cm_rf2)

# 2 stressmodels
ml3.get_step_response("river").plot()
ml3.get_step_response("river2").plot()

plt.legend(
[
"Before (ChangeModel)",
"After (ChangeModel)",
"Before (method 2)",
"After (method 2)",
]
)

[9]:

<matplotlib.legend.Legend at 0x7fbfbc37d690>


# References#

Obergfell, C., Bakker, M. and Maas, K. (2019), Identification and Explanation of a Change in the Groundwater Regime using Time Series Analysis. Groundwater, 57: 886-894. https://doi.org/10.1111/gwat.12891