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, three recharge models are currently available:

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.


To run this notebook and the related non-linear recharge models, it is strongly recommended to install Numba ( This Just-In-Time (JIT) compiler compiles the computationally intensive part of the recharge calculation, making the non-linear model as fast as the Linear recharge model.

import pandas as pd
import pastas as ps
import matplotlib.pyplot as plt

Python version: 3.7.8 | packaged by conda-forge | (default, Jul 31 2020, 02:37:09)
[Clang 10.0.1 ]
Numpy version: 1.18.5
Scipy version: 1.4.0
Pandas version: 1.1.2
Pastas version: 0.16.0b
Matplotlib version: 3.1.3
numba version: 0.49.1

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=True)

# Make this millimeters per day
evap = ps.read_knmi("../data/etmgeg_260.txt", variables="EV24").series * 1e3
rain = ps.read_knmi("../data/etmgeg_260.txt", variables="RH").series * 1e3

fig, axes = plt.subplots(3,1, figsize=(10,6), sharex=True)
head.plot(ax=axes[0], x_compat=True, linestyle=" ", marker=".")
evap.plot(ax=axes[1], x_compat=True)
rain.plot(ax=axes[2], x_compat=True)
axes[0].set_ylabel("Head [m]")
axes[1].set_ylabel("Evap [mm/d]")
axes[2].set_ylabel("Rain [mm/d]")

plt.xlim("1985", "2005");
INFO: Inferred frequency for time series EV24 260: freq=D
INFO: Inferred frequency for time series RH 260: freq=D

Make a basic model

The normal workflow may be used to create and calibrate the model. 1. Create a Pastas Model instance 2. Choose a recharge model. All recharge models can be accessed through the recharge subpackage (ps.rch). 3. Create a RechargeModel object and add it to the model 4. Solve and visualize the model

ml = ps.Model(head)

# Select a recharge model
rch = ps.rch.FlexModel()
#rch = ps.rch.Berendrecht()
#rch = ps.rch.Linear()

rm = ps.RechargeModel(rain, evap, recharge=rch, rfunc=ps.Gamma, name="rch")

ml.solve(noise=True, tmin="1990", report="basic")
INFO: Cannot determine frequency of series head: freq=None. The time series is irregular.
INFO: Inferred frequency for time series RH 260: freq=D
INFO: Inferred frequency for time series EV24 260: freq=D
Fit report head                   Fit Statistics
nfev     27                     EVP           89.46
nobs     351                    R2             0.89
noise    True                   RMSE           0.07
tmin     1990-01-01 00:00:00    AIC           21.06
tmax     2005-10-14 00:00:00    BIC           63.53
freq     D                      Obj            0.47
warmup   3650 days 00:00:00     ___
solver   LeastSquares           Interpolated     No

Parameters (8 were optimized)
                optimal   stderr     initial   vary
rch_A          0.427248   ±6.00%    0.721210   True
rch_n          0.668178   ±2.96%    1.000000   True
rch_a        300.022264  ±13.50%   10.000000   True
rch_srmax     54.702385   ±5.59%  250.000000   True
rch_lp         0.250000    ±nan%    0.250000  False
rch_ks        19.547379  ±11.15%  100.000000   True
rch_gamma      4.051118   ±9.43%    4.000000   True
rch_simax      2.000000    ±nan%    2.000000  False
rch_kv         1.000000    ±nan%    1.000000  False
constant_d     0.805622   ±4.06%    1.359779   True
noise_alpha   34.051304  ±11.71%   15.000000   True

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("A").sum()
ax =,3))
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 (ml.solve(noise=False)) and then solve it again using a noise model.

  • For relatively shallow groundwater levels, it may be better to use the Exponential response function as the the non-linear models already cause a delayed response.


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 - The precipitation and evapotranspiration (etmgeg_260.txt) are downloaded from