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 numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pastas as ps

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

/tmp/ipykernel_948/1496237964.py:2: DeprecationWarning: 
Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),
(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)
but was not found to be installed on your system.
If this would cause problems for you,
please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466
        
  import pandas as pd

Python version: 3.11.6
NumPy version: 1.26.4
Pandas version: 2.2.0
SciPy version: 1.12.0
Matplotlib version: 3.8.3
Numba version: 0.59.0
LMfit version: 1.2.2
Latexify version: Not Installed
Pastas version: 1.4.0

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

../_images/1b4a53a45d641b10fea9cedc0aa385f3c442011d694034cc6185d5aa8dd4f2e9.png

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 0x7f9acf107310>

../_images/4b06c9e0d16c5cc697a3d96d9464263a737f7297d01a009caa26144184673b67.png

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

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", noise=True)
ml.plots.results(figsize=(10, 6))

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

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", noise=True)
ml2.plots.results(figsize=(10, 6));

/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/timeseries_utils.py:90: FutureWarning: Day.delta is deprecated and will be removed in a future version. Use pd.Timedelta(obj) instead
  if hasattr(offset, "delta"):

/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/timeseries_utils.py:90: FutureWarning: Day.delta is deprecated and will be removed in a future version. Use pd.Timedelta(obj) instead
  if hasattr(offset, "delta"):

../_images/3cc5d3b03bc07acdbba02b9460c38de35e49cf3b88bfd9114ed36319adb2b057.png

../_images/17a6b376d54545ac1110cf331d540b4d5abc8ce4be6c54249ab8dd2534ff2b6e.png

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	dist	stderr	optimal
test_1_A	1.00000	1.000000e-05	100.0	True	test	uniform	0.014107	0.449464
test_1_a	10.00000	1.000000e-02	1000.0	True	test	uniform	0.137745	2.466723
test_2_A	1.00000	1.000000e-05	100.0	True	test	uniform	0.009413	0.851243
test_2_a	10.00000	1.000000e-02	1000.0	True	test	uniform	0.045022	1.124797
test_beta	0.00000	-inf	inf	True	test	uniform	6.555924	-49.877036
test_tchange	734503.00000	6.124110e+05	825914.0	False	test	uniform	NaN	734503.000000
step_d	1.00000	NaN	NaN	True	step	uniform	0.019171	1.029063
step_tstart	734503.00000	6.124110e+05	825914.0	False	step	uniform	NaN	734503.000000
constant_d	330.18797	NaN	NaN	True	constant	uniform	0.011379	329.301851
noise_alpha	7.00000	1.000000e-05	5000.0	True	noise	uniform	1.727050	28.265550

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"]);

../_images/e4fa8f1f173baff5d933da77d6a9967a1c43b5c558e5fffa24d6ee199d9d27fe.png

5. Another way#

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

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

/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/timeseries_utils.py:90: FutureWarning: Day.delta is deprecated and will be removed in a future version. Use pd.Timedelta(obj) instead
  if hasattr(offset, "delta"):
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/timeseries_utils.py:90: FutureWarning: Day.delta is deprecated and will be removed in a future version. Use pd.Timedelta(obj) instead
  if hasattr(offset, "delta"):

../_images/d81535783a1ef85e9df4fd36d5d6a2b9172a1cead584fcce50cdbbc3866f8b43.png

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 0x7f9acd0a8890>

../_images/ad0a1d14f55c0e75b5da9c7e85b4462f45b02fad6f7625387ab8930584ffb50d.png

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