Example 2: Analysis of groundwater monitoring networks using Pastas#

This notebook is supplementary material to the following paper in 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.6
Numpy      : 2.5.2
Pandas     : 3.0.5
Scipy      : 1.18.1
Matplotlib : 3.11.1
Numba      : 0.67.0

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

../../../_images/d01f0739fb18eb7718236c856e283f96cd2964ee396843eb568af757a898c6f2.png

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
    ps.RechargeModel(ml, rain, evap, rfunc=ps.Gamma(), name="recharge")

    # Add a StressModel to simulate the effect of the groundwater extractions
    sm = ps.StressModel(
        ml, well / 1e6, rfunc=ps.Hantush(), name="well", settings="well", up=False
    )

    # Add a NoiseModel (explicitly required since Pastas 1.5)
    nm = ps.ArNoiseModel(ml)

    # Estimate the model parameters
    ps.solver.Lmfit(ml)
    ml.solve(tmin=tmin, tmax=tmax, report=False)

    # 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 29
   nm = ps.ArNoiseModel(ml)

   # Estimate the model parameters
   ps.solver.Lmfit(ml)
---> 29     ml.solve(tmin=tmin, tmax=tmax, report=False)

   # Check if the estimated effect of the groundwater extraction is significant.
   # If not, delete the stressmodel and calibrate the model again.

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/model.py:1132, in Model.solve(self, tmin, tmax, freq, warmup, solver, report, initial, weights, fit_constant, freq_obs, initialize, reset_settings, noise, **kwargs)
   LeastSquares(model=self)
# Solve model
-> 1132 solve_success, result = self.solver.solve(weights=weights, **kwargs)
# Update the parameters with the results from the optimization
for column in result.columns:

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/solver/least_squares.py:1166, in Lmfit.solve(self, weights, **kwargs)
objfunction = partial(
   self.objfunction,
   noise=noise,
   weights=weights,
)
mini = lmfit.Minimizer(
   userfcn=objfunction,
   calc_covar=True,
   params=parameters,
   **kwargs,
)
-> 1166 self.result = mini.minimize(method=self.method)
names = self.result.var_names
# Set all parameter attributes

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/lmfit/minimizer.py:2355, in Minimizer.minimize(self, method, params, **kws)
       if (key.lower().startswith(user_method) or
               val.lower().startswith(user_method)):
           kwargs['method'] = val
-> 2355 return function(**kwargs)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/lmfit/minimizer.py:1674, in Minimizer.leastsq(self, params, max_nfev, **kws)
result.call_kws = lskws
try:
-> 1674     lsout = scipy_leastsq(self.__residual, variables, **lskws)
except AbortFitException:
   pass

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/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)
   if maxfev == 0:
       maxfev = 200*(n + 1)
--> 439     retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
                            gtol, maxfev, epsfcn, factor, diag)
else:
   if col_deriv:

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/lmfit/minimizer.py:540, in Minimizer.__residual(self, fvars, apply_bounds_transformation)
   self.result.success = False
   raise AbortFitException(f"fit aborted: too many function evaluations {self.max_nfev}")
--> 540 out = self.userfcn(params, *self.userargs, **self.userkws)
if callable(self.iter_cb):
   abort = self.iter_cb(params, self.result.nfev, out,
                        *self.userargs, **self.userkws)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/solver/least_squares.py:1209, in Lmfit.objfunction(self, parameters, noise, weights)
"""Objective function that is minimized by the Lmfit solver."""
p = np.array([p.value for p in parameters.values()])
-> 1209 return misfit(
   model=self.model,
   p=p,
   noise=noise,
   weights=weights,
   callback=None,
   returnseparate=False,
)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/solver/objective/misfit.py:45, in misfit(model, p, noise, weights, callback, returnseparate)
"""
Shared objective function for solvers to calculate residuals or noise.

   (...)     42     The calculated residuals or noise, optionally with separate components.
"""
# Get the residuals or the noise
---> 45 res = model.residuals(p)
if noise:
   res = model.noise(p=p, res=res) * model._noise_weights(p=p, res=res)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/model.py:723, in Model.residuals(self, p, tmin, tmax, freq, warmup)
"""Calculate the residual series.

Parameters
   (...)    720     Series with the residuals.
"""
obs = self.observations(tmin=tmin, tmax=tmax, freq=freq)
--> 723 sim = self._simulate_on_observations(
   p=p, tmin=tmin, tmax=tmax, freq=freq, warmup=warmup
)
res = obs.subtract(sim).rename("Residuals")
if res.hasnans:

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/model.py:635, in Model._simulate_on_observations(self, p, tmin, tmax, freq, warmup)
freq_obs = freq if settings["freq_obs"] is None else settings["freq_obs"]
# simulate model
--> 635 sim = self.simulate(
   p=p, tmin=tmin, tmax=tmax, freq=freq, warmup=warmup, return_warmup=False
)
# Get the oseries calibration series
obs = self.observations(tmin=tmin, tmax=tmax, freq=freq_obs)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/model.py:551, in Model.simulate(self, p, tmin, tmax, freq, warmup, return_warmup)
for sm in self.stressmodels.values():
   p_sm = p[istart : istart + sm.nparam]
--> 551     contrib = sm.simulate(
       p=p_sm,
       tmin=sim_index_min,
       tmax=tmax,
       freq=freq,
       dt=dt,
   )
   if contrib.hasnans:
       logger.error(
           f"StressModel {sm.name} contribution simulation"
           f" with parameters {p_sm} contains NaN-values."
       )

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/stressmodels.py:611, in StressModel.simulate(self, p, tmin, tmax, freq, dt)
def simulate(
   self,
   p: ArrayLike,
   (...)    608     dt: float = 1.0,
) -> Series:
   """Simulate the stressmodel's contribution."""
--> 611     return self._simulate(tuple(p), tmin, tmax, freq, dt)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/decorators.py:391, in conditional_cachedmethod.<locals>.decorator.<locals>.wrapper(self, *args, **kwargs)
   return cached_func(self, *args, **kwargs)
else:
--> 391     return func(self, *args, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/pastas/stressmodels.py:654, in StressModel._simulate(self, p, tmin, tmax, freq, dt)
stress = self.stress.series
npoints = stress.index.size
h = Series(
--> 654     data=fftconvolve(stress, b, "full")[:npoints],
   index=stress.index,
   name=self.name,
   dtype=np.asarray(p).dtype,
)
return h

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/scipy/signal/_signaltools.py:712, in fftconvolve(in1, in2, mode, axes)
s2 = in2.shape
shape = [max((s1[i], s2[i])) if i not in axes else s1[i] + s2[i] - 1
        for i in range(in1.ndim)]
--> 712 ret = _freq_domain_conv(xp, in1, in2, axes, shape, calc_fast_len=True)
return _apply_conv_mode(ret, s1, s2, mode, axes, xp=xp)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/scipy/signal/_signaltools.py:540, in _freq_domain_conv(xp, in1, in2, axes, shape, calc_fast_len)
   in2 = xp.astype(in2, xp_default_dtype(xp))
sp1 = fft(in1, fshape, axes=axes)
--> 540 sp2 = fft(in2, fshape, axes=axes)
ret = ifft(sp1 * sp2, fshape, axes=axes)
if calc_fast_len:

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/scipy/fft/_backend.py:29, in _ScipyBackend.__ua_function__(method, args, kwargs)
if fn is None:
   return NotImplemented
---> 29 return fn(*args, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/scipy/fft/_basic_backend.py:138, in rfftn(x, s, axes, norm, overwrite_x, workers, plan)
def rfftn(x, s=None, axes=None, norm=None,
         overwrite_x=False, workers=None, *, plan=None):
--> 138     return _execute_nD('rfftn', _duccfft.rfftn, x, s=s, axes=axes, norm=norm,
                      overwrite_x=overwrite_x, workers=workers, plan=plan)

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/scipy/fft/_basic_backend.py:57, in _execute_nD(func_str, duccfft_func, x, s, axes, norm, overwrite_x, workers, plan)
if is_numpy(xp):
   x = np.asarray(x)
---> 57     return duccfft_func(x, s=s, axes=axes, norm=norm,
                         overwrite_x=overwrite_x, workers=workers, plan=plan)
norm = _validate_fft_args(workers, plan, norm)
if hasattr(xp, 'fft'):

File ~/checkouts/readthedocs.org/user_builds/pastas/envs/stable/lib/python3.14/site-packages/scipy/fft/_duccfft/basic.py:177, in r2cn(forward, x, s, axes, norm, overwrite_x, workers, plan)
   raise ValueError("at least 1 axis must be transformed")
# Note: overwrite_x is not utilized
--> 177 return pfft.r2c(tmp, axes, forward, norm, None, workers)

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=list(mls.keys()), columns=["optimal", "stderr"], dtype=float
)
for name, ml in mls.items():
    if "well" in ml.stressmodels:
        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)

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

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