Study area and data

As an example to illustrate the methodology, this study focuses on four locations with a groundwater monitoring well in the northern part of the Dutch province Limburg, see Fig 2. The wells are located in the management area of regional water authority Waterschap Limburg. These wells have limited external influences, such as limited interaction with surface water or major abstractions. The northern part of Limburg has sandy soils which are susceptible to groundwater droughts due to their reliance on precipitation for groundwater recharge (Brakkee et al., Reference Brakkee, van Huijgevoort and Bartholomeus2021; van den Eertwegh et al., Reference van den Eertwegh, de Louw, Witte, van Huijgevoort, Bartholomeus, van Deijl, van Dam, Hunink, America, Pouwels, Hoefsloot and de Wit2021).

Figure 2.

The four monitoring wells and the KNMI meteorological stations in northern Limburg (KNMI, 2023b, 2023a). The monitoring wells are located near Ell, Heibloem, Sevenum and Mariapeel.

Three data types are used in this study:

• Observed phreatic groundwater level time series. These time series are used as the goal function to train the TFN model. The four monitoring wells provide phreatic groundwater observations for the period 2012–2020;

• Observed meteorological time series (precipitation and reference crop evapotranspiration). This time series is used as input to train the TFN model (von Asmuth et al., Reference von Asmuth, Baggelaar, Bakker, Brakenhoff, Collenteur, Ebbens, Mondeel, Klop and Schaars2021). The observed daily precipitation sum is obtained from nearby KNMI weather stations Ell, Heibloem and Sevenum. KNMI station Ell provides daily observed reference crop evapotranspiration data. In the remainder of this article, we refer to the reference crop evapotranspiration as evapotranspiration data;

• Detrended long-term meteorological time series (precipitation and evapotranspiration). These time series are provided for the period 1910–2022 and are detrended for climatological trends based on the last 30 years. They represent the meteorological conditions over a long time period in the current climate. These time series are generated by performing a seasonal detrending approach using a LOESS-fit on long-term homogeneous observations from various KNMI meteorological stations in the Netherlands. This approach accounts for the average seasonal trend change. The detrended series are then gridded by performing a spatial interpolation: ordinary kriging for the precipitation data and thin plate spline for the evapotranspiration data. More details can be found in (STOWA, 2023). The climatological properties of historical meteorological data are validated. The Kolmogorov–Smirnov test revealed that the historical precipitation at locations Heibloem and Sevenum revealed deviations from the current climate. Therefore, the cumulative density functions of these precipitation series were bias corrected using quantile mapping. As a result, the historical meteorological data does match the current climate and can be used to simulate groundwater levels in the current climate.

Table 1 shows an overview of the data and their characteristics. The locations of the point observations are visualised in Fig 2.

Table 1.

The type, temporal resolution, and time period of the data

Transfer function-noise modelling

We set up a TFN-model using the open-source Python package Pastas (Collenteur et al., Reference Collenteur, Bakker, Caljé, Klop and Schaars2019). Pastas utilises the Predefined Impulse Response Functions in Continuous Time (PIRFICT) methodology (von Asmuth et al., Reference von Asmuth, Bierkens and Maas2002). Effectively, the change in groundwater levels is described by

(1) $${h\left( t \right) = \mathop \sum \limits_{m = 1}^M {h_m}\left( t \right) + d + r\left( t \right)\;}$$

In which h(t) is the observed groundwater level [m] at time t [day], h _m (t) is the contribution of input stress to the groundwater level [m] at time t, d is the base level of the model [m], and r(t) is the model residual [m] at time t (Collenteur et al., Reference Collenteur, Bakker, Caljé, Klop and Schaars2019).

The selected locations have limited external influences, according to Waterschap Limburg. Hence, it is assumed that the main explanatory variables for groundwater dynamics are precipitation and evapotranspiration, ignoring external factors such as surface water level dynamics and groundwater abstractions. We define daily groundwater recharge as input stress, which consists of the daily precipitation minus daily evapotranspiration. The contribution of the recharge input stress to the groundwater level can be determined using the following impulse response function:

(2) $${{h_{recharge}}\left( t \right) = \mathop \int \nolimits_{ - \infty }^t {S_{recharge}}\left( \tau \right){\vartheta _{recharge}}\left( {t - \tau } \right)d\tau }$$

In which S _recharge is the time series of the recharge input stress at preceding time τ [m] and ϑ _m is the impulse response function of the recharge input stress [−] (von Asmuth et al., Reference von Asmuth, Bierkens and Maas2002). The impulse response function describes the behaviour of groundwater levels following an impulse of the input stress. A commonly used impulse response function is a function based on the gamma distribution (Besbes & De Marsily, Reference Besbes and De Marsily1984), which has parameters for the shape (a [day], n [−]) and scaling factor (A [day⁻ⁿ]):

(3) $${\vartheta \left( t \right) = {{A{t^{n - 1}}{e^{ - {t \over a}}}\;} \over {{a^n}{\rm{\Gamma }}\left( n \right)}}}$$

Another important aspect of time series modelling is the noise model, which is used to reduce autocorrelation in the model residuals and simulate stochastic system behaviour. The residual difference between observed and simulated groundwater levels reflects system behaviour that cannot be described by the input stress(es). Model residuals often show a correlation that makes them dependent on the previous day’s residual. When this correlation effect persists for several days, it potentially leads to exaggerated simulated groundwater levels. To address this issue, correlated noise is introduced using an exponential decay (Collenteur et al., Reference Collenteur, Bakker, Caljé, Klop and Schaars2019):

(4) $$n_{c}\left(t\right)=e^{-{\Delta t \over \alpha }}\cdot n_{c}\left(t-1\right)+n\left(t\right)$$

In which n _c is correlated noise [−], α is the noise decay parameter [day] and n is uncorrelated (white) noise result of a random process [−]. The created noise series, with a time step Δt [day] of one day, is by definition the difference between the observations and the deterministic part (Eq. 1).

Pastas uses both observed groundwater time series and the observed meteorological series to calibrate the parameters of the transfer function, using the non-linear least squares method for the optimisation. The model is trained for an eight year-period. This time period is chosen such that the model represents the current response of the groundwater system (Zaadnoordijk et al., Reference Zaadnoordijk, Bus, Lourens and Berendrecht2019), assuming that no significant alterations have taken place in the water system. The calibration period is six years (2014−2020) and the validation period is two years (2012−2014). The calibration period consists of both dry and moderate wet years to include these dynamics in the calibration procedure. The validation describes model performance for a different period than the calibration period.

The model performance is evaluated using two goodness of fit metrics: the explained variance percentage (EVP) and the root mean square error (RMSE). These metrics are commonly used for groundwater level TFN models (von Asmuth et al., Reference von Asmuth, Maas, Knotters, Bierkens, Bakker, Olsthoorn, Cirkel, Leunk, Schaars and von Asmuth2012). The EVP [%] describes the percentage of variation in the groundwater levels that is explained by precipitation and potential evapotranspiration:

(5) $${EVP = {{{\sigma _h}^2 - {\sigma _r}^2} \over {{\sigma _h}^2}}*100\%}$$

where σ_h ² [m²] is the variance of the groundwater observations and σ_r ² [m²] is the variance of the residuals. The returned value is bounded between 0 and 100%, higher percentages indicate more explained variance and therefore a better fit. The RMSE [m] measures the average difference between simulated and observed values and is calculated by the square root of the mean of the residuals:

(6) $$RMSE=\sqrt{{\sum r^{2} \over N_{r}}}$$

where N _r [−] is the amount of residuals r [m]. A lower root mean square error implies smaller residuals and therefore a better fit. The goodness of fit criteria for the model performance are set to a minimum EVP of 70% and RMSE below 0.5, which often serve as a rule of thumb for TFN models (Pezij et al., Reference Pezij, Augustijn, Hendriks, Hulscher, S., Kew, van Der Wiel, Wanders and van Oldenborgh2020; van Engelenburg et al., Reference van Engelenburg, de Jonge, Rijpkema, van Slobbe and Bense2020; Ritter et al, Reference Ritter and Muñoz-Carpena2013).

Characterising groundwater droughts

Many methods to describe groundwater droughts exist, most notably using percentiles (Zaadnoordijk et al., Reference Zaadnoordijk, Bus, Lourens and Berendrecht2019) or the Standardised Groundwater Index (Bloomfield & Marchant, Reference Bloomfield and Marchant2013). In this paper, we use annual minimum groundwater levels to characterise the intensity and frequency of groundwater droughts. This indicator is chosen for its simplicity to illustrate the approach, it is not necessarily the most appropriate indicator to describe groundwater droughts. Therefore, the duration of drought is not considered. The drought intensity is represented by the absolute value of the annual minimum groundwater level [m +NAP], while the frequency of droughts is assessed using return periods for the annual minima [years], estimated using their plotting position (Benard & Bos-Levenbach, Reference Benard and Bos-Levenbach1954):

(7) $$F_{i}={i-0.3 \over N+0.4}$$

In which F_i is the plotting position for a specific drought intensity on a frequency curve [year⁻¹], i denotes the rank number of the drought intensity [−], sorted from high to low annual minimum water levels [m +NAP], where the lowest value has a rank of 1. N is the total number of data values [−]. The return period T [year] is the inverse of the survival function of F_i . The return period indicates how often a drought event of a specific intensity can be expected to occur:

(8) $$T={1 \over 1-F_{i}}$$

The intensity and frequency of groundwater droughts are derived from the long-term generated groundwater levels (1910−2022) and compared to the intensity and frequency of groundwater droughts derived from observations of the last eight years.