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.6
Numpy : 2.4.6
Pandas : 3.0.3
Scipy : 1.18.0
Matplotlib : 3.11.0
Numba : 0.66.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");
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
ps.solver.Lmfit(model=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}")
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of RechargeModel needs to be a Pastas Model. Please provide the model as the first argument: RechargeModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of StressModel needs to be a Pastas Model. Please provide the model as the first argument: StressModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of ArNoiseModel needs to be a Pastas Model. Please provide the model as the first argument: ArNoiseModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_noisemodel is deprecated and will not be available from Pastas version >= 2.4.0. Noise models are now added by adding the Pastas Model as the first argument during noise model initialization (i.e., ps.ArNoiseModel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of RechargeModel needs to be a Pastas Model. Please provide the model as the first argument: RechargeModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of StressModel needs to be a Pastas Model. Please provide the model as the first argument: StressModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of ArNoiseModel needs to be a Pastas Model. Please provide the model as the first argument: ArNoiseModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_noisemodel is deprecated and will not be available from Pastas version >= 2.4.0. Noise models are now added by adding the Pastas Model as the first argument during noise model initialization (i.e., ps.ArNoiseModel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of RechargeModel needs to be a Pastas Model. Please provide the model as the first argument: RechargeModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of StressModel needs to be a Pastas Model. Please provide the model as the first argument: StressModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of ArNoiseModel needs to be a Pastas Model. Please provide the model as the first argument: ArNoiseModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_noisemodel is deprecated and will not be available from Pastas version >= 2.4.0. Noise models are now added by adding the Pastas Model as the first argument during noise model initialization (i.e., ps.ArNoiseModel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of RechargeModel needs to be a Pastas Model. Please provide the model as the first argument: RechargeModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of StressModel needs to be a Pastas Model. Please provide the model as the first argument: StressModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_stressmodel is deprecated and will not be available from Pastas version >= 2.4.0. Stressmodels are now added by adding the Pastas Model as the first argument during stressmodel initialization (i.e., ps.Stressmodel(model=ml, *args))
warn(message=msg, category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:140: FutureWarning: From Pastas 2.4, the first argument of ArNoiseModel needs to be a Pastas Model. Please provide the model as the first argument: ArNoiseModel(model=ml, ...).
warn(message=msg % (self._name, self._name), category=FutureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:212: FutureWarning: add_noisemodel is deprecated and will not be available from Pastas version >= 2.4.0. Noise models are now added by adding the Pastas Model as the first argument during noise model initialization (i.e., ps.ArNoiseModel(model=ml, *args))
warn(message=msg, category=FutureWarning)
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[4], line 32
28 ml.add_noisemodel(nm)
29
30 # Estimate the model parameters
31 ps.solver.Lmfit(model=ml)
---> 32 ml.solve(tmin=tmin, tmax=tmax, report=False)
33
34 # Check if the estimated effect of the groundwater extraction is significant.
35 # If not, delete the stressmodel and calibrate the model again.
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:1017, in Model.solve(self, tmin, tmax, freq, warmup, solver, report, initial, weights, fit_constant, freq_obs, initialize, reset_settings, noise, **kwargs)
1014 LeastSquares(model=self)
1016 # Solve model
-> 1017 solve_success, result = self.solver.solve(weights=weights, **kwargs)
1018 # Update the parameters with the results from the optimization
1019 for column in result.columns:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/solver/least_squares.py:1159, in Lmfit.solve(self, noise, weights, **kwargs)
1148 objfunction = partial(
1149 self.objfunction,
1150 noise=noise,
1151 weights=weights,
1152 )
1153 mini = lmfit.Minimizer(
1154 userfcn=objfunction,
1155 calc_covar=True,
1156 params=parameters,
1157 **kwargs,
1158 )
-> 1159 self.result = mini.minimize(method=self.method)
1160 names = self.result.var_names
1162 # Set all parameter attributes
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/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/latest/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/latest/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/latest/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/latest/lib/python3.14/site-packages/pastas/solver/least_squares.py:1202, in Lmfit.objfunction(self, parameters, noise, weights)
1200 """Objective function that is minimized by the Lmfit solver."""
1201 p = np.array([p.value for p in parameters.values()])
-> 1202 return misfit(
1203 ml=self.model,
1204 p=p,
1205 noise=noise,
1206 weights=weights,
1207 callback=None,
1208 returnseparate=False,
1209 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/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/latest/lib/python3.14/site-packages/pastas/model.py:727, in Model._noise_weights(self, p, tmin, tmax, freq, warmup)
724 p = self.get_parameters()
726 # Calculate the residuals
--> 727 res = self.residuals(p, tmin, tmax, freq, warmup)
729 # Calculate the weights
730 weights = self.noisemodel.weights(res, p[-self.noisemodel.nparam :])
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:611, in Model.residuals(self, p, tmin, tmax, freq, warmup)
606 freq_obs = (
607 freq if self.settings["freq_obs"] is None else self.settings["freq_obs"]
608 )
610 # simulate model
--> 611 sim = self.simulate(
612 p=p, tmin=tmin, tmax=tmax, freq=freq, warmup=warmup, return_warmup=False
613 )
615 # Get the oseries calibration series
616 obs = self.observations(tmin=tmin, tmax=tmax, freq=freq_obs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:531, in Model.simulate(self, p, tmin, tmax, freq, warmup, return_warmup)
529 istart = 0 # Track parameters index to pass to stressmodel object
530 for sm in self.stressmodels.values():
--> 531 contrib = sm.simulate(
532 p=p[istart : istart + sm.nparam],
533 tmin=sim_index.min(),
534 tmax=tmax,
535 freq=freq,
536 dt=dt,
537 )
538 sim = sim.add(contrib)
539 istart += sm.nparam
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:612, in StressModel.simulate(self, p, tmin, tmax, freq, dt)
603 def simulate(
604 self,
605 p: ArrayLike,
(...) 609 dt: float = 1.0,
610 ) -> Series:
611 """Simulate the stressmodel's contribution."""
--> 612 return self._simulate(tuple(p), tmin, tmax, freq, dt)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/decorators.py:342, in conditional_cachedmethod.<locals>.decorator.<locals>.wrapper(self, *args, **kwargs)
340 return cached_func.__get__(self, type(self))(*args, **kwargs)
341 else:
--> 342 return func(self, *args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:651, in StressModel._simulate(self, p, tmin, tmax, freq, dt)
623 """Simulate the head contribution.
624
625 Parameters
(...) 648
649 """
650 self.update_stress(tmin=tmin, tmax=tmax, freq=freq)
--> 651 b = self._get_block(p, dt, tmin, tmax)
652 stress = self.stress.series
653 npoints = stress.index.size
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/stressmodels.py:350, in StressModelBase._get_block(self, p, dt, tmin, tmax)
348 else:
349 maxtmax = None
--> 350 b = self.rfunc.block(p, dt, maxtmax=maxtmax)
351 return b
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/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/latest/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/latest/lib/python3.14/site-packages/pastas/rfunc.py:1274, in Hantush.step(self, p, dt, cutoff, maxtmax, **kwargs)
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/latest/lib/python3.14/site-packages/pastas/rfunc.py:1191, in Hantush.numpy_step(A, a, b, t)
1189 tau2 = tau[inv_mask]
1190 b_tau2 = b_over_tau[inv_mask]
-> 1191 F[inv_mask] = 2.0 * k0rho - w * exp1(tau2) + w_minus_1 * exp1(tau2 + b_tau2)
1193 return A * F / (2.0 * k0rho)
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))
