Example 2: Analysis of groundwater monitoring networks using Pastas#
This notebook is supplementary material to the following paper submitted to Groundwater:
Collenteur, R.A., Bakker, M., Caljé, R., Klop, S.A., Schaars, F. (2019) Pastas: open source software for the analysis of groundwater time series. Groundwater. doi: 10.1111/gwat.12925.
In this second example, it is demonstrated how scripts can be used to analyze a large number of time series. Consider a pumping well field surrounded by a number of observations wells. The pumping wells are screened in the middle aquifer of a three-aquifer system. The objective is to estimate the drawdown caused by the groundwater pumping in each observation well.
1. Import the packages#
# Import the packages
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pastas as ps
ps.show_versions()
ps.set_log_level("ERROR")
try:
from timml import ModelMaq, Well
plot_timml = True
except ImportError:
plot_timml = False
plot_results = False
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
2. Importing the time series#
In this codeblock the time series are imported. The following time series are imported:
44 time series with head observations [m] from the monitoring network;
precipitation [m/d] from KNMI station Oudenbosch;
potential evaporation [m/d] from KNMI station de Bilt;
Total pumping rate [m3/d] from well field Seppe.
# Dictionary to hold all heads
heads = {}
# Load a metadata-file with xy-coordinates from the groundwater heads
metadata_heads = pd.read_csv("data/metadata_heads.csv", index_col=0)
distances = pd.read_csv("data/distances.csv", index_col=0)
# Add the groundwater head observations to the database
for fname in os.listdir("./data/heads/"):
fname = os.path.join("./data/heads/", fname)
obs = pd.read_csv(fname, parse_dates=True, index_col=0).squeeze()
heads[obs.name] = obs
# Load a metadata-file with xy-coordinates from the explanatory variables
metadata = pd.read_csv("data/metadata_stresses.csv", index_col=0)
# Import the precipitation, evaporation and well time series
rain = pd.read_csv("data/rain.csv", parse_dates=True, index_col=0).squeeze()
evap = pd.read_csv("data/evap.csv", parse_dates=True, index_col=0).squeeze()
well = pd.read_csv("data/well.csv", parse_dates=True, index_col=0).squeeze()
# Plot the stresses
fig, [ax1, ax2, ax3] = plt.subplots(3, 1, figsize=(10, 5), sharex=True)
rain.plot(ax=ax1)
evap.plot(ax=ax2)
well.plot(ax=ax3)
plt.xlim("1960", "2018");
3/4/5. Creating and optimizing the Time Series Model#
For each time series of groundwater head observations a TFN model is constructed with the following model components:
A Constant
A NoiseModel
A RechargeModel object to simulate the effect of recharge
A StressModel object to simulate the effect of groundwater extraction
Calibrating all models can take a couple of minutes!!
# Create folder to save the model figures
mls = {}
mlpath = "models"
if not os.path.exists(mlpath):
os.mkdir(mlpath)
# Choose the calibration period
tmin = "1970"
tmax = "2017-09"
num = 0
for name, head in heads.items():
# Create a Model for each time series and add a StressModel2 for the recharge
ml = ps.Model(head, name=name)
# Add the RechargeModel to simulate the effect of rainfall and evaporation
rm = ps.RechargeModel(rain, evap, rfunc=ps.Gamma(), name="recharge")
ml.add_stressmodel(rm)
# Add a StressModel to simulate the effect of the groundwater extractions
sm = ps.StressModel(
well / 1e6, rfunc=ps.Hantush(), name="well", settings="well", up=False
)
ml.add_stressmodel(sm)
# Add a NoiseModel (explicitly required since Pastas 1.5)
nm = ps.ArNoiseModel()
ml.add_noisemodel(nm)
# Estimate the model parameters
ml.solve(tmin=tmin, tmax=tmax, report=False, solver=ps.solver.Lmfit())
# Check if the estimated effect of the groundwater extraction is significant.
# If not, delete the stressmodel and calibrate the model again.
gain, stderr = ml.parameters.loc["well_A", ["optimal", "stderr"]]
if stderr is None:
stderr = 10.0
if 1.96 * stderr > -gain:
num += 1
ml.del_stressmodel("well")
ml.solve(tmin=tmin, tmax=tmax, report=False)
# Plot the results and store the plot
mls[name] = ml
if plot_results:
ml.plots.results()
path = os.path.join(mlpath, name + ".png")
plt.savefig(path, bbox_inches="tight")
plt.close()
print(f"The number of models where the well is dropped from the model is {num}")
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[4], line 31
27 nm = ps.ArNoiseModel()
28 ml.add_noisemodel(nm)
29
30 # Estimate the model parameters
---> 31 ml.solve(tmin=tmin, tmax=tmax, report=False, solver=ps.solver.Lmfit())
32
33 # Check if the estimated effect of the groundwater extraction is significant.
34 # If not, delete the stressmodel and calibrate the model again.
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/model.py:990, in Model.solve(self, tmin, tmax, freq, warmup, solver, report, initial, weights, fit_constant, freq_obs, initialize, reset_settings, noise, **kwargs)
987 self.add_solver(solver=LeastSquares())
989 # Solve model
--> 990 solve_success, result = self.solver.solve(weights=weights, **kwargs)
991 # Update the parameters with the results from the optimization
992 for column in result.columns:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/solver/least_squares.py:1157, in Lmfit.solve(self, noise, weights, **kwargs)
1146 objfunction = partial(
1147 self.objfunction,
1148 noise=noise,
1149 weights=weights,
1150 )
1151 mini = lmfit.Minimizer(
1152 userfcn=objfunction,
1153 calc_covar=True,
1154 params=parameters,
1155 **kwargs,
1156 )
-> 1157 self.result = mini.minimize(method=self.method)
1158 names = self.result.var_names
1160 # Set all parameter attributes
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/lmfit/minimizer.py:2355, in Minimizer.minimize(self, method, params, **kws)
2352 if (key.lower().startswith(user_method) or
2353 val.lower().startswith(user_method)):
2354 kwargs['method'] = val
-> 2355 return function(**kwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/lmfit/minimizer.py:1674, in Minimizer.leastsq(self, params, max_nfev, **kws)
1672 result.call_kws = lskws
1673 try:
-> 1674 lsout = scipy_leastsq(self.__residual, variables, **lskws)
1675 except AbortFitException:
1676 pass
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/scipy/optimize/_minpack_py.py:439, in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
437 if maxfev == 0:
438 maxfev = 200*(n + 1)
--> 439 retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
440 gtol, maxfev, epsfcn, factor, diag)
441 else:
442 if col_deriv:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/lmfit/minimizer.py:540, in Minimizer.__residual(self, fvars, apply_bounds_transformation)
537 self.result.success = False
538 raise AbortFitException(f"fit aborted: too many function evaluations {self.max_nfev}")
--> 540 out = self.userfcn(params, *self.userargs, **self.userkws)
542 if callable(self.iter_cb):
543 abort = self.iter_cb(params, self.result.nfev, out,
544 *self.userargs, **self.userkws)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/solver/least_squares.py:1200, in Lmfit.objfunction(self, parameters, noise, weights)
1198 """Objective function that is minimized by the Lmfit solver."""
1199 p = np.array([p.value for p in parameters.values()])
-> 1200 return misfit(
1201 ml=self.ml,
1202 p=p,
1203 noise=noise,
1204 weights=weights,
1205 callback=None,
1206 returnseparate=False,
1207 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/solver/objective_function.py:43, in misfit(ml, p, noise, weights, callback, returnseparate)
41 # Get the residuals or the noise
42 if noise:
---> 43 rv = ml.noise(p) * ml._noise_weights(p)
44 else:
45 rv = ml.residuals(p)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/model.py:687, in Model.noise(self, p, tmin, tmax, freq, warmup)
684 p = self.get_parameters()
686 # Calculate the residuals
--> 687 res = self.residuals(p, tmin, tmax, freq, warmup)
688 p = p[-self.noisemodel.nparam :]
690 # Calculate the noise
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/model.py:592, in Model.residuals(self, p, tmin, tmax, freq, warmup)
587 freq_obs = (
588 freq if self.settings["freq_obs"] is None else self.settings["freq_obs"]
589 )
591 # simulate model
--> 592 sim = self.simulate(
593 p=p, tmin=tmin, tmax=tmax, freq=freq, warmup=warmup, return_warmup=False
594 )
596 # Get the oseries calibration series
597 obs = self.observations(tmin=tmin, tmax=tmax, freq=freq_obs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/model.py:512, in Model.simulate(self, p, tmin, tmax, freq, warmup, return_warmup)
510 istart = 0 # Track parameters index to pass to stressmodel object
511 for sm in self.stressmodels.values():
--> 512 contrib = sm.simulate(
513 p=p[istart : istart + sm.nparam],
514 tmin=sim_index.min(),
515 tmax=tmax,
516 freq=freq,
517 dt=dt,
518 )
519 sim = sim.add(contrib)
520 istart += sm.nparam
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/stressmodels.py:565, in StressModel.simulate(self, p, tmin, tmax, freq, dt)
556 def simulate(
557 self,
558 p: ArrayLike,
(...) 562 dt: float = 1.0,
563 ) -> Series:
564 """Simulate the stressmodel's contribution."""
--> 565 return self._simulate(tuple(p), tmin, tmax, freq, dt)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/decorators.py:324, in conditional_cachedmethod.<locals>.decorator.<locals>.wrapper(self, *args, **kwargs)
322 return cached_func.__get__(self, type(self))(*args, **kwargs)
323 else:
--> 324 return func(self, *args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/stressmodels.py:604, in StressModel._simulate(self, p, tmin, tmax, freq, dt)
576 """Simulate the head contribution.
577
578 Parameters
(...) 601
602 """
603 self.update_stress(tmin=tmin, tmax=tmax, freq=freq)
--> 604 b = self._get_block(p, dt, tmin, tmax)
605 stress = self.stress.series
606 npoints = stress.index.size
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/stressmodels.py:344, in StressModelBase._get_block(self, p, dt, tmin, tmax)
342 else:
343 maxtmax = None
--> 344 b = self.rfunc.block(p, dt, maxtmax=maxtmax)
345 return b
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/rfunc.py:358, in RfuncBase.block(self, p, dt, cutoff, maxtmax, **kwargs)
334 """Return the block function.
335
336 Parameters
(...) 355 Block response values.
356 """
357 if self.use_block:
--> 358 b = self.block_from_step(
359 p=p, dt=dt, cutoff=cutoff, maxtmax=maxtmax, **kwargs
360 )
361 else:
362 b = self.block_from_impulse(
363 p=p, dt=dt, cutoff=cutoff, maxtmax=maxtmax, **kwargs
364 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/rfunc.py:430, in RfuncBase.block_from_step(self, p, dt, cutoff, maxtmax, **kwargs)
400 def block_from_step(
401 self,
402 p: ArrayLike,
(...) 406 **kwargs,
407 ) -> ArrayLike:
408 """Return the block function from the step response.
409
410 Parameters
(...) 428 Block response values.
429 """
--> 430 s = self.step(p=p, dt=dt, cutoff=cutoff, maxtmax=maxtmax, **kwargs)
431 b = np.append(s[0], np.subtract(s[1:], s[:-1]))
432 return b
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/rfunc.py:1269, in Hantush.step(self, p, dt, cutoff, maxtmax, **kwargs)
1248 """Return the step function for Hantush response.
1249
1250 Parameters
(...) 1266 Array with the step response.
1267 """
1268 A, a, b = p
-> 1269 t = self.get_t(p=p, dt=dt, cutoff=cutoff, maxtmax=maxtmax, **kwargs)
1271 step = (
1272 self.quad_step(A=A, a=a, b=b, t=t)
1273 if self.quad
1274 else self.numpy_step(A=A, a=a, b=b, t=t)
1275 )
1276 return step
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/dev/lib/python3.14/site-packages/pastas/rfunc.py:316, in RfuncBase.get_t(self, p, dt, cutoff, maxtmax, **kwargs)
284 def get_t(
285 self,
286 p: ArrayLike,
(...) 290 **kwargs,
291 ) -> ArrayLike:
292 """Determine the times at which to evaluate the step response, from t=0.
293
294 Parameters
(...) 314 Times at which the response is evaluated.
315 """
--> 316 if np.ndim(dt) > 0:
317 return np.asarray(dt, dtype=float)
319 tmax = self._resolve_tmax(p=np.asarray(p).real, cutoff=cutoff, **kwargs)
KeyboardInterrupt:
Make plots for publication#
In the next codeblocks the Figures used in the Pastas paper are created. The following figures are created:
Figure of the drawdown estimated for each observations well;
Figure of the decomposition of the different contributions;
Figure of the pumping rate of the well field.
Figure of the drawdown estimated for each observations well#
x = np.linspace(100, 5000, 100)
if plot_timml:
# Values from REGIS II v2.2 (Site id B49F0240)
z = [9, -25, -83, -115, -190] # Reference to NAP
kv = np.array(
[
1e-3,
5e-3,
]
) # Min-Max of Vertical hydraulic conductivity for both leaky layer
D1 = z[0] - z[1] # Estimated thickness of leaky layer
c1 = D1 / kv # Estimated resistance
D2 = z[2] - z[3]
c2 = D2 / kv
kh1 = np.array(
[
1e0,
2.5e0,
]
) # Min-Max of Horizontal hydraulic conductivity for aquifer 1
kh2 = np.array(
[
1e1,
2.5e1,
]
) # Min-Max of Horizontal hydraulic conductivity for aquifer 2
mlm = ModelMaq(
kaq=[kh1.mean(), 35], z=z, c=[c1.max(), c2.mean()], topboundary="semi", hstar=0
)
w = Well(mlm, 0, 0, 34791, layers=1)
mlm.solve()
h = mlm.headalongline(x, 0)
np.savetxt("head_timml.out", h)
else:
h = np.loadtxt("head_timml.out")
# Get the parameters and distances to plot
params = pd.DataFrame(index=mls.keys(), columns=["optimal", "stderr"], dtype=float)
for name, ml in mls.items():
if "well" in ml.stressmodels.keys():
params.loc[name] = (
ml.parameters.loc["well_A", ["optimal", "stderr"]]
* well.loc["2007":].mean()
/ 1e6
)
# Select model per aquifer
shallow = metadata_heads.z.loc[(metadata_heads.z < 96)].index
aquifer = metadata_heads.z.loc[(metadata_heads.z < 186) & (metadata_heads.z > 96)].index
# Make the plot
fig = plt.figure(figsize=(8, 5))
plt.grid(zorder=-10)
display_error_bars = True
if display_error_bars:
plt.errorbar(
distances.loc[shallow, "Seppe"],
params.loc[shallow, "optimal"],
yerr=1.96 * params.loc[shallow, "stderr"],
linestyle="",
elinewidth=2,
marker="",
markersize=10,
capsize=4,
)
plt.errorbar(
distances.loc[aquifer, "Seppe"],
params.loc[aquifer, "optimal"],
yerr=1.96 * params.loc[aquifer, "stderr"],
linestyle="",
elinewidth=2,
marker="",
capsize=4,
)
plt.scatter(
distances.loc[shallow],
params.loc[shallow, "optimal"],
marker="^",
s=80,
label="aquifer 1",
)
plt.scatter(
distances.loc[aquifer],
params.loc[aquifer, "optimal"],
marker="s",
s=80,
label="aquifer 2",
)
# Plot two-layer TimML model for comparison
plt.plot(x, h[0], color="C0", linestyle="--", label="TimML L1")
plt.plot(x, h[1], color="C1", linestyle="--", label="TimML L2")
plt.ylabel("steady drawdown (m)")
plt.xlabel("radial distance from the center of the well field (m)")
plt.xlim(0, 4501)
plt.legend(loc=4)
<matplotlib.legend.Legend at 0x7d8456930980>
Example figure of a TFN model#
# Select a model to plot
ml = mls["B49F0232_5"]
# Create the figure
[ax1, ax2, ax3] = ml.plots.decomposition(
split_contributions=False, figsize=(7, 6), ytick_base=1, tmin="1985"
)
plt.xticks(rotation=0)
ax1.set_yticks([2, 0, -2])
ax1.set_ylabel("head (m)")
ax1.legend().set_visible(False)
ax3.set_yticks([-4, -6])
ax2.set_ylabel(
"contributions (m) "
) # Little trick to get the label right
ax3.set_xlabel("year")
ax3.set_ylabel("")
ax3.set_title("pumping well")
Text(0.5, 1.0, 'pumping well')
Figure of the pumping rate of the well field#
fig, ax = plt.subplots(1, 1, figsize=(8, 2.5), sharex=True)
ax.plot(well, color="k")
ax.set_ylabel("pumping rate\n[m$^3$/day]")
ax.set_xlabel("year")
ax.set_xlim(pd.Timestamp("1951"), pd.Timestamp("2018"))
(np.float64(-6940.0), np.float64(17532.0))
