Non-linear recharge models#
R.A. Collenteur, University of Graz
This notebook explains the use of the RechargeModel stress model to simulate the combined effect of precipitation and potential evaporation on the groundwater levels. For the computation of the groundwater recharge, four recharge models are currently available:
Berendrecht(Berendrecht et al., 2006)FlexModel(Collenteur et al., 2021)Peterson(Peterson & Western, 2014)
The first model is a simple linear function of precipitation and potential evaporation while the latter two are simulate a non-linear response of recharge to precipitation using a soil-water balance concepts. Detailed descriptions of these models can be found in articles listed in the References at the end of this notebook.
import matplotlib.pyplot as plt
import pandas as pd
import pastas as ps
ps.show_versions()
ps.set_log_level("INFO")
Pastas : 2.0.0
Python : 3.14.6
Numpy : 2.4.6
Pandas : 3.0.5
Scipy : 1.18.0
Matplotlib : 3.11.1
Numba : 0.66.0
Read Input data#
Input data handling is similar to other stressmodels. The only thing that is necessary to check is that the precipitation and evaporation are provided in mm/day. This is necessary because the parameters for the non-linear recharge models are defined in mm for the length unit and days for the time unit. It is possible to use other units, but this would require manually setting the initial values and parameter boundaries for the recharge models.
head = pd.read_csv(
"data/B32C0639001.csv", parse_dates=["date"], index_col="date"
).squeeze()
# Make this millimeters per day
evap = pd.read_csv("data/evap_260.csv", index_col=0, parse_dates=[0]).squeeze()
rain = pd.read_csv("data/rain_260.csv", index_col=0, parse_dates=[0]).squeeze()
ps.plots.series(
head,
[evap, rain],
figsize=(10, 6),
labels=["Head [m]", "Evap [mm/d]", "Rain [mm/d]"],
);
Make a basic model#
The normal workflow may be used to create and calibrate the model.
Create a Pastas
ModelinstanceChoose a recharge model. All recharge models can be accessed through the recharge subpackage (
ps.rch).Create a
RechargeModelobject and add it to the modelSolve and visualize the model
ml = ps.Model(head)
ps.ArNoiseModel(model=ml)
# Select a recharge model
rch = ps.rch.FlexModel(gw_uptake=True)
# rch = ps.rch.Berendrecht()
# rch = ps.rch.Linear()
# rch = ps.rch.Peterson()
rm = ps.RechargeModel(ml, rain, evap, recharge=rch, rfunc=ps.Gamma(), name="rch")
ml.solve(tmin="1990", report="basic")
ml.plots.results(figsize=(10, 6));
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[3], line 12
8 # rch = ps.rch.Peterson()
9
10 rm = ps.RechargeModel(ml, rain, evap, recharge=rch, rfunc=ps.Gamma(), name="rch")
11
---> 12 ml.solve(tmin="1990", report="basic")
13 ml.plots.results(figsize=(10, 6));
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:830, in LeastSquares.solve(self, weights, **kwargs)
815 bounds = Bounds(
816 lb=pmin,
817 ub=pmax,
818 keep_feasible=True,
819 )
821 objfunction = partial(
822 self.objfunction,
823 noise=noise,
(...) 827 callback=kwargs.pop("callback", None),
828 )
--> 830 self.result = least_squares(
831 fun=objfunction,
832 x0=initial[vary],
833 jac=self.jac,
834 bounds=bounds,
835 method=self.method,
836 ftol=self.ftol,
837 xtol=self.xtol,
838 gtol=self.gtol,
839 x_scale=self.x_scale,
840 loss=self.loss,
841 f_scale=self.f_scale,
842 max_nfev=self.max_nfev,
843 diff_step=self.diff_step,
844 tr_solver=self.tr_solver,
845 **kwargs,
846 )
848 self.pcov = DataFrame(
849 self.get_covariances(
850 self.result.jac,
(...) 856 columns=parameters.index,
857 )
859 # Prepare return values
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/_lib/_util.py:660, in _workers_wrapper.<locals>.inner(*args, **kwds)
658 with MapWrapper(_workers) as mf:
659 kwargs['workers'] = mf
--> 660 return func(*args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/least_squares.py:1019, in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs, callback, workers)
1015 result = call_minpack(vector_fun.fun, x0, vector_fun.jac, ftol, xtol, gtol,
1016 max_nfev, x_scale, jac_method=jac)
1018 elif method == 'trf':
-> 1019 result = trf(vector_fun.fun, vector_fun.jac, x0, f0, J0, lb, ub, ftol, xtol,
1020 gtol, max_nfev, x_scale, loss_function, tr_solver,
1021 tr_options.copy(), verbose, callback=callback_wrapped)
1023 elif method == 'dogbox':
1024 if tr_solver == 'lsmr' and 'regularize' in tr_options:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/trf.py:124, in trf(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, callback)
120 return trf_no_bounds(
121 fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev, x_scale,
122 loss_function, tr_solver, tr_options, verbose, callback=callback)
123 else:
--> 124 return trf_bounds(
125 fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
126 loss_function, tr_solver, tr_options, verbose, callback=callback)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/trf.py:376, in trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, callback)
372 f_true = f.copy()
374 cost = cost_new
--> 376 J = jac(x)
377 njev += 1
379 if loss_function is not None:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:750, in VectorFunction.jac(self, x)
748 def jac(self, x):
749 self._update_x(x)
--> 750 self._update_jac()
751 if hasattr(self.J, "astype"):
752 # returns a copy so that downstream can't overwrite the
753 # internal attribute. But one can't copy a LinearOperator
754 return self.J.astype(self.J.dtype)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:719, in VectorFunction._update_jac(self)
716 else:
717 self._njev += 1
--> 719 self.J = self.jac_wrapped(xp_copy(self.x), f0=self.f)
720 self.J_updated = True
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:454, in _VectorJacWrapper.__call__(self, x, f0, **kwds)
452 self.njev += 1
453 elif self.jac in FD_METHODS:
--> 454 J, dct = approx_derivative(
455 self.fun,
456 x,
457 f0=f0,
458 **self.finite_diff_options,
459 )
460 self.nfev += dct['nfev']
462 if self.sparse_jacobian:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_numdiff.py:611, in approx_derivative(fun, x0, method, rel_step, abs_step, f0, bounds, sparsity, as_linear_operator, args, kwargs, full_output, workers)
609 with MapWrapper(workers) as mf:
610 if sparsity is None:
--> 611 J, _nfev = _dense_difference(fun_wrapped, x0, f0, h,
612 use_one_sided, method,
613 mf)
614 else:
615 if not issparse(sparsity) and len(sparsity) == 2:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_numdiff.py:738, in _dense_difference(fun, x0, f0, h, use_one_sided, method, workers)
735 l = next(gen)
736 u = next(gen)
--> 738 f1 = next(f_evals)
739 f2 = next(f_evals)
740 if one_sided:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_numdiff.py:912, in _Fun_Wrapper.__call__(self, x)
909 if xp.isdtype(x.dtype, "real floating"):
910 x = xp.astype(x, self.x0.dtype)
--> 912 f = np.atleast_1d(self.fun(x, *self.args, **self.kwargs))
913 if f.ndim > 1:
914 raise RuntimeError("`fun` return value has "
915 "more than 1 dimension.")
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:424, in _VectorFunWrapper.__call__(self, x)
422 def __call__(self, x):
423 self.nfev += 1
--> 424 return np.atleast_1d(self.fun(x))
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/least_squares.py:263, in _WrapArgsKwargs.__call__(self, x)
262 def __call__(self, x):
--> 263 return self.f(x, *self.args, **self.kwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/solver/least_squares.py:772, in LeastSquares.objfunction(self, p, noise, weights, initial, vary, callback)
770 par = initial
771 par[vary] = p
--> 772 return misfit(
773 ml=self.model, p=par, noise=noise, weights=weights, callback=callback
774 )
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:706, in Model.noise(self, p, tmin, tmax, freq, warmup)
703 p = self.get_parameters()
705 # Calculate the residuals
--> 706 res = self.residuals(p, tmin, tmax, freq, warmup)
707 p = p[-self.noisemodel.nparam :]
709 # Calculate the noise
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:509, in Model.simulate(self, p, tmin, tmax, freq, warmup, return_warmup)
501 # Get the simulation index and the time step
502 # Check if the requested index matches the model settings
503 if (
504 tmin == self.settings["tmin"]
505 and tmax == self.settings["tmax"]
506 and freq == self.settings["freq"]
507 and warmup == self.settings["warmup"]
508 ):
--> 509 sim_index = self.sim_index
510 else:
511 # simulate with the requested settings, but do not update
512 # the model settings, since this is just for one time
513 sim_index = _get_sim_index(
514 tmin=tmin - warmup,
515 tmax=tmax,
516 freq=freq,
517 time_offset=self.time_offset,
518 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:1389, in Model.sim_index(self)
1371 @property
1372 def sim_index(self) -> DatetimeIndex:
1373 """Property that returns the simulation index, including the warmup.
1374
1375 Using the tmin, tmax, freq, and warmup from the model
(...) 1383 model is simulated.
1384 """
1385 return _get_sim_index(
1386 tmin=self.settings["tmin"] - self.settings["warmup"],
1387 tmax=self.settings["tmax"],
1388 freq=self.settings["freq"],
-> 1389 time_offset=self.time_offset,
1390 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:1358, in Model.time_offset(self)
1355 if st.freq_original:
1356 # calculate the offset from the default frequency
1357 t = st.series_original.index
-> 1358 base = t.min().ceil(freq)
1359 mask = t >= base
1360 if np.any(mask):
File pandas/_libs/tslibs/timestamps.pyx:3082, in pandas._libs.tslibs.timestamps.Timestamp.ceil()
-> 3082 'Could not get source, probably due dynamically evaluated source code.'
File pandas/_libs/tslibs/timestamps.pyx:2754, in pandas._libs.tslibs.timestamps.Timestamp._round()
-> 2754 'Could not get source, probably due dynamically evaluated source code.'
File pandas/_libs/tslibs/offsets.pyx:6362, in pandas._libs.tslibs.offsets.to_offset()
-> 6362 'Could not get source, probably due dynamically evaluated source code.'
File pandas/_libs/tslibs/offsets.pyx:1288, in pandas._libs.tslibs.offsets.Tick.__mul__()
-> 1288 'Could not get source, probably due dynamically evaluated source code.'
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/numpy/_core/numeric.py:2494, in isclose(a, b, rtol, atol, equal_nan)
2491 y = float(y)
2493 # atol and rtol can be arrays
-> 2494 if not (np.all(np.isfinite(atol)) and np.all(np.isfinite(rtol))):
2495 err_s = np.geterr()["invalid"]
2496 err_msg = f"One of rtol or atol is not valid, atol: {atol}, rtol: {rtol}"
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/numpy/_core/fromnumeric.py:2634, in all(a, axis, out, keepdims, where)
2548 @array_function_dispatch(_all_dispatcher)
2549 def all(a, axis=None, out=None, keepdims=np._NoValue, *, where=np._NoValue):
2550 """
2551 Test whether all array elements along a given axis evaluate to True.
2552
(...) 2632
2633 """
-> 2634 return _wrapreduction_any_all(a, np.logical_and, 'all', axis, out,
2635 keepdims=keepdims, where=where)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/numpy/_core/fromnumeric.py:97, in _wrapreduction_any_all(obj, ufunc, method, axis, out, **kwargs)
95 pass
96 else:
---> 97 return reduction(axis=axis, out=out, **passkwargs)
99 return ufunc.reduce(obj, axis, bool, out, **passkwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/numpy/_core/_methods.py:70, in _all(a, axis, dtype, out, keepdims, where)
68 # Parsing keyword arguments is currently fairly slow, so avoid it for now
69 if where is True:
---> 70 return umr_all(a, axis, dtype, out, keepdims)
71 return umr_all(a, axis, dtype, out, keepdims, where=where)
KeyboardInterrupt:
Analyze the estimated recharge flux#
After the parameter estimation we can take a look at the recharge flux computed by the model. The flux is easy to obtain using the get_stress method of the model object, which automatically provides the optimal parameter values that were just estimated. After this, we can for example look at the yearly recharge flux estimated by the Pastas model.
recharge = ml.get_stress("rch").resample("YE").sum()
ax = recharge.plot.bar(figsize=(10, 3))
ax.set_xticklabels(recharge.index.year)
plt.ylabel("Recharge [mm/year]");
A few things to keep in mind:#
Below are a few things to keep in mind while using the (non-linear) recharge models.
The use of an appropriate warmup period is necessary, so make sure the precipitation and evaporation are available some time (e.g., one year) before the calibration period.
Make sure that the units of the precipitation fluxes are in mm/day and that the DatetimeIndex matches exactly.
It may be possible to fix or vary certain parameters, dependent on the problem. Obtaining better initial parameters may be possible by solving without a noise model first and then solve it again using a noise model.
For relatively shallow groundwater levels (e.g., < 10 meters), it may be better to use the
Exponentialresponse function as the the non-linear models already cause a delayed response.
References#
Berendrecht, W. L., Heemink, A. W., van Geer, F. C., and Gehrels, J. C. (2003), Decoupling of modeling and measuring interval in groundwater time series analysis based on response characteristics, Journal of Hydrology, 278, 1–16.
Berendrecht, W. L., Heemink, A. W., van Geer, F. C., and Gehrels, J. C. (2006), A non-linear state space approach to model groundwater fluctuations, Advances in Water Resources, 29, 959–973.
Collenteur, R. A., Bakker, M., Klammler, G., and Birk, S. (2021), Estimation of groundwater recharge from groundwater levels using nonlinear transfer function noise models and comparison to lysimeter data, Hydrol. Earth Syst. Sci., 25, 2931–2949.
Peterson, T.J., and A.W. Western (2014), Nonlinear time-series modeling of unconfined groundwater head, Water Resour. Res., 50, 8330-8355.
Von Asmuth, J.R., Maas, K., Bakker, M. and Petersen, J. (2008), Modeling Time Series of Ground Water Head Fluctuations Subjected to Multiple Stresses. Groundwater, 46: 30-40.
Data Sources#
In this notebook we analysed a head time series near the town of De Bilt in the Netherlands. Data is obtained from the following resources:
The heads (
B32C0639001.csv) are downloaded from https://www.dinoloket.nl/The precipitation and evapotranspiration (
rain_260.csvandevap_260.csv) are downloaded from https://knmi.nl