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.

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

import pastas as ps

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

Pastas     : 2.0.0
Python     : 3.14.0
Numpy      : 2.4.6
Pandas     : 3.0.3
Scipy      : 1.17.1
Matplotlib : 3.10.9
Numba      : 0.65.1

1. Load the data

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

axes = ps.plots.series(
    head=head,
    stresses=[prec, evap, river],
    tmin="2000",
    labels=["Head\n[m]", "Prec\n[mm/d]", "Evap\n[mm/d]", "River [m3/d]"],
)

2. The weighting factor

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

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()

<matplotlib.legend.Legend at 0x7a882931f380>

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.

# Normal Model
ml = ps.Model(head, name="linear")
ml.add_noisemodel(ps.ArNoiseModel())

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

ml.add_stressmodel([sm, step])
ml.solve(report=False, tmin="2004", tmax="2017-12-31")
ml.plots.results(figsize=(10, 6))

# ChangeModel
ml2 = ps.Model(head, name="linear")
ml2.add_noisemodel(ps.ArNoiseModel())

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

ml2.add_stressmodel([cm, step])
ml2.solve(report=False, tmin="2004", tmax="2017-12-31")
ml2.plots.results(figsize=(10, 6));

/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: RuntimeWarning: overflow encountered in exp
  omega = 1 / (np.exp(beta * (t - sigma)) + 1)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:2704: 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.

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.14

ml2.parameters

	initial	pmin	pmax	vary	name	optimal	stderr
test_1_A	1.00000	1.000000e-05	100.0	True	test	0.851244	0.020368
test_1_a	10.00000	1.000000e-02	10000.0	True	test	1.124806	0.097426
test_2_A	1.00000	1.000000e-05	100.0	True	test	0.449465	0.030528
test_2_a	10.00000	1.000000e-02	10000.0	True	test	2.466751	0.298081
test_beta	0.00000	-inf	inf	True	test	49.875242	14.186233
test_tchange	734503.00000	6.124110e+05	825914.0	False	test	734503.000000	NaN
step_A	1.00000	NaN	NaN	True	step	1.029062	0.041485
step_tstart	734503.00000	6.124110e+05	825914.0	False	step	734503.000000	NaN
constant_d	330.18797	NaN	NaN	True	constant	329.301851	0.024624
noise_alpha	7.00000	1.000000e-05	5000.0	True	noise	28.265548	3.737270

4. Compare the response functions

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

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.

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

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.add_stressmodel([r1, r2, step])
ml3.solve(report=False, tmin="2004", tmax="2017-12-31")
ml3.plots.results(figsize=(10, 6));

How do the results compare?

# 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)",
    ]
)

<matplotlib.legend.Legend at 0x7a88200e0050>

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