Introduction
The availability of groundwater is important for many purposes in the Netherlands, such as agriculture, industry, drinking water, prevention of land subsidence and mitigation of salinisation problems. However, during periods of precipitation deficit (meteorological drought), groundwater levels can drop significantly, resulting in groundwater droughts and reduced water availability for many of these user functions (van Loon, Reference van Loon2015). Recent droughts have received much attention due to the severe damage they have caused to a wide range of sectors (Brakkee et al., Reference Brakkee, van Huijgevoort and Bartholomeus2021). For example, the drought of 2018 caused significant damage to the agricultural sector, hampered navigation, drinking water supplies and disrupted industry, with estimated economic losses ranging from €450 to €2080 million (van Hussen et al., Reference van Hussen, van De Velde, Läkamp and van Der Kooij2019).
Several factors, including a changing climate, increased groundwater abstractions and changes in the water system, are putting pressure on the groundwater availability. The increase in climate variability is expected to lead to even drier summers in the future (Philip et al., Reference Philip, Kew, van Der Wiel, Wanders and van Oldenborgh2020). The recently published KNMI’23 climate projections confirm this trend, predicting drier summers in the Netherlands under all expected future climate scenarios (van Dorland et al., Reference van Dorland, Beersma, Bessembinder, Bloemendaal, van Den Brink, Brotons Blanes, Drijfhout, Groenland, Haarsma, Homan, Keizer, Krikken, Le Bars, Lenderink, van Meijgaard, Meirink, Overbeek, Reerink, Selten and van Der Wiel2024). Population growth, economic growth and changes in land use increase the demand for groundwater and the corresponding increase in groundwater abstractions in dry years (CBS et al., 2024). Moreover, the current design and use of the Dutch water system have left areas susceptible to droughts. The accelerated discharge of precipitation in wet winter periods through drainage pipes, ditches and (channelised) streams plays a major role in this (de Lenne & Worm, Reference de Lenne and Worm2020).
Monitoring groundwater dynamics is essential to gain insight into groundwater droughts. Groundwater observations can serve as the basis for calculating drought indicators, and various statistical methods are used to characterise current groundwater conditions using these observations. Characterisation of groundwater droughts involves providing insight into their characteristics such as intensity, duration and frequency (Petersen-Perlman et al., Reference Petersen-Perlman, Aguilar-Barajas and Megdal2022). Ideally, at least 30 years of observations are required to capture long-term groundwater trends (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). Detection of extreme droughts requires an even longer period. However, the length of groundwater observations useful for analysis is often limited to 5–10 years. Reliable long-term groundwater measurements (over a period of more than 50 years) are hardly available in the Netherlands (Verhagen & Avis, Reference Verhagen and Avis2021). In addition, these historical groundwater observations are often not applicable for long-term analyses due to changes in climate, abstractions and water systems (Ritzema et al., Reference Ritzema, Heuvelink, Heinen, Bogaart, van der Bolt, Hack-ten Broeke, Hoogland, Knotters, Massop and Vroon2012; 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). Therefore, it is challenging to accurately characterise groundwater droughts due to this lack of long-term groundwater observations. Extended time series are essential to improve the characterisation of groundwater droughts.
Groundwater modelling can be used in combination with meteorological data to extend groundwater time series. Numerical models, such as MODFLOW, can be used for this purpose. However, these models are computationally expensive, which limits their use for long-term modelling. An alternative approach is time series modelling, which is commonly used to simulate groundwater levels (Bakker & Schaars, Reference Bakker and Schaars2019).
This study applies a transfer function-noise (TFN) model, which uses regression analysis to establish relationships between observed groundwater levels and various stress factors, such as precipitation, evaporation, surface water levels and groundwater abstractions. This relationship can be used to extend or interpolate observed time series of groundwater levels. This data-driven approach is often simpler and much faster than applying a numerical groundwater model (Bakker & Schaars, Reference Bakker and Schaars2019). Furthermore, the response of a water system to input stresses is included, which improves our understanding of the effects of an input stress (Pezij et al., Reference Pezij, Augustijn, Hendriks, Hulscher, S., Kew, van Der Wiel, Wanders and van Oldenborgh2020). Finally, the stochastic system dynamics are modelled using a noise model (von Asmuth et al., Reference von Asmuth, Bierkens and Maas2002). Many studies have shown the effectiveness of using TFN modelling in groundwater studies (Collenteur et al., Reference Collenteur, Bakker, Klammler and Birk2021, Reference Collenteur, Moeck, Schirmer and Birk2023; Jemeljanova et al., Reference Jemeljanova, Collenteur, Kmoch, Bikše, Popovs and Kalvans2023; Pezij et al., Reference Pezij, Augustijn, Hendriks, Hulscher, S., Kew, van Der Wiel, Wanders and van Oldenborgh2020; Rudolph et al., Reference Rudolph, Collenteur, Kavousi, Giese, Wöhling, Birk, Hartmann and Reimann2023; Zaadnoordijk et al., Reference Zaadnoordijk, Bus, Lourens and Berendrecht2019).
The extension of observed groundwater data using precipitation and evaporation is constrained by the availability of long-term meteorological measurements reflective of the current climate. Due to climate variability, only 30 years of meteorological data can be used to extend groundwater time series, as this period defines the current climate (World Meteorological Organization, 2017). Therefore, the approach needs to move towards climate models or artificially generated meteorological data for further extension. Furthermore, the use of synthetic meteorological data in time series modelling shows potential for extending groundwater level time series (El Mezouary et al., Reference El Mezouary, El Mansouri and El Bouhaddioui2020; Vonk, Reference Vonk2021). This study uses detrended historical meteorological time series (STOWA, 2023).
Therefore, we propose a methodology to generate long-term groundwater level time series by combining groundwater modelling and historical detrended meteorological time series for the period 1910–2022. The aim of this methodology is to improve the understanding of the intensity and occurrence of groundwater droughts.
Methods
Figure 1 outlines the used methodology. First, a TFN model is developed using both observed groundwater levels and observed meteorological data. The trained model is then used to generate long-term groundwater levels for the period 1910–2022 using the detrended historical meteorological data. The intensity and frequency of groundwater droughts are then characterised using annual minimum groundwater levels for both the generated long-term groundwater simulations and short-term groundwater measurements. Finally, a comparison is made between groundwater droughts occurring in the simulated long-term groundwater levels and the observations. Each of these stages is described in detail in the following sections of this chapter.
Figure 1. Flowchart of methodology.
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).
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
$${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:
$${{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−n]):
$${\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):
$$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:
$${EVP = {{{\sigma _h}^2 - {\sigma _r}^2} \over {{\sigma _h}^2}}*100\%}$$
where σh 2 [m2] is the variance of the groundwater observations and σr 2 [m2] 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:
$$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):
$$F_{i}={i-0.3 \over N+0.4}$$
In which Fi is the plotting position for a specific drought intensity on a frequency curve [year−1], 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 Fi . The return period indicates how often a drought event of a specific intensity can be expected to occur:
$$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.
Results
Calibration and validation of the groundwater simulations
First, we present the results of the model calibrations and validations to ensure that the long-term generated groundwater levels are representative of the current climate and groundwater system. Figure 3 shows the calibration and validation results for the four locations. The TFN model is able to capture seasonal groundwater dynamics. In addition, the model is able to capture extreme events, as visible during the dry summer period of 2018. The goodness of fit metrics are shown in the figures. These metrics meet the criteria for the calibration period as well as the validation period. The calibrated model parameters can be found in the Appendix.
Figure 3. Calibration (orange) and validation (blue) of the time series model, together with goodness of fit metrics EVP and RMSE with respect to measurements (dots).
Long-term groundwater simulation (1910–2022)
We then used the trained TFN-models in combination with the long-term detrended historical meteorological time series to generate long-term groundwater simulations. Figure 4 shows the resulting groundwater levels for the period 1910−2022. It must be emphasised that these results represent an estimate of what groundwater levels would have been under the current climate and water system, rather than reflecting actual observed historical groundwater levels. The results are historical projections based on current conditions, so they are useful for describing current groundwater conditions.
Figure 4. Long-term groundwater levels based on simulations with the calibrated model and detrended historical meteorological data. The solid red line is the recent dry year 2018, which can be compared to observed groundwater levels. The dashed red line is the record-year 1976, which is the year with the largest observed precipitation deficit.
As can be seen, groundwater levels in 1976 are very low, even lower than in 2018. In fact, there are several years that were drier than 2018, a trend that is consistent across all locations.
A validation process was carried out to ensure the accuracy of the long-term groundwater level simulations. First, long-term generated groundwater simulations were compared with actual groundwater observations within the analysis period (2012−2020). Similar performance on goodness-of-fit metrics confirmed the model’s ability to reproduce observed groundwater levels for the calibration period, as shown in Fig 5.
Figure 5. Validation of long-term generated groundwater levels (green) using measurements in the period 2012–2020. The figure also includes the results of the model calibration (orange) and validation (blue).
Comparing groundwater droughts
We then used annual minimum groundwater levels to characterise groundwater droughts for the long-term generated groundwater simulations (1910−2022) and groundwater observations of the last eight years. Figure 6 shows the intensity of groundwater levels, characterised by annual minimum levels, and the frequency, represented by the return periods of these annual minimum levels.
Figure 6. Differences in return periods for annual minimum groundwater levels.
The use of long-term groundwater simulations allows the characterisation of historical droughts, such as the one of 1976. In particular, the long-term simulations show that several years, namely 1921 and 1976, experienced a greater drought intensity than 2018, when the groundwater levels are adjusted to the current climate. Thus, based on the generated long-term groundwater levels, the minimum value in 2018 appears comparatively less severe than based on observations from the last eight years.
Analysing the intensity and frequency of groundwater droughts using both the long-term generated groundwater levels and the observations from the last eight years reveals substantial differences. Figure 6 shows a significant contrast in the frequency of groundwater droughts between long-term simulations and observations. This difference is evident in the different distribution functions of groundwater droughts shown in the figure, implying that long-term simulations yield different return periods compared to groundwater observations.
Table 2 shows the results using annual minimum groundwater levels to characterise the 2018 drought. While the intensity is consistent across observations and simulations, the results for frequency differ for the groundwater observations and the long-term groundwater simulations. Specifically, the frequency derived from the simulations indicates a return period that is larger than the return period derived using observations.
Table 2. Difference in the characterisation of the annual minimum water level for reference year 2018 between groundwater observations and simulations
Discussion
The characterisation of groundwater droughts reveals substantial differences in return periods derived from groundwater observations and long-term simulations. The methodology presented in this study allows a direct conversion of groundwater levels into return periods and vice versa. Consequently, the use of long-term generated groundwater levels for the current climate provides more statistically accurate estimates of drought frequency (Fig 7a). Case in point is the 2018 drought, which occurs less frequently than expected based on observations (once in 12 years) compared to the simulations (once in 24 years). The change in the frequency distribution also affects the groundwater level corresponding to a given return period (Fig 7b). For example, the groundwater levels corresponding to a drought that occurs once every 10 years (T10) is much lower when derived from observations (26.48 m + NAP) compared to derived from the simulations (26.78 m +NAP). Arguably, the method thus facilitates a more accurate quantification of past drought intensity under current climate and groundwater system.
Figure 7. Differences in return periods between groundwater observations (2012–2020) and long-term simulations (1910–2022).
The results reveal a significant change in the frequency patterns of long-term simulated groundwater levels near the T10 return period. Prior to this change, groundwater levels generally decline gradually with increasing return periods, followed by a sudden drop before stabilising. This trend is observed at all locations, although for some locations this effect is more pronounced. Possible reasons for this include physical processes that are noticeable at very low groundwater levels, or inaccuracies in the estimation of return times. Estimation of return periods assumes independence between events, affecting the statistical analysis. As a result, the estimation of return periods is influenced by the selection of drought events, which can lead to bias. In addition, longer return periods are less accurate because they have occurred less in the simulations. Therefore, with 112 years of generated data, reliable return periods can be estimated up to about T30, which occur on average more than four times in the long-term generated groundwater levels. Additionally, some groundwater levels are briefly simulated above the surface level. This anomaly is likely due to the model being calibrated during a relatively dry period (2014−2020), resulting in less accurate modelling of periods with high precipitation. Lastly, the serial correlation of the residuals in the model should be considered. If the model’s residuals are not randomly distributed, it can lead to systematic deviations over long periods, affecting the model’s reliability over extended timescales (Zaadnoordijk, Reference Zaadnoordijk2022).
There are also practical ramifications of using our proposed approach. Water management in the Netherlands often follows a risk-based approach, taking into account the frequency of extreme events. The differences in the results of this study have significant implications, suggesting that reliance on observation-based criteria alone may lead to underestimation of drought intensity and consequently delay necessary action. In addition, other indicators or statistics can be used to assess the intensity and frequency of low groundwater levels, which can help to define criteria for the implementation of measures such as irrigation restrictions.
A limitation of the study is that characterising groundwater levels using annual minimum groundwater levels and estimating return periods using plotting positions is a shallow approach and not always indicative of drought conditions. However, the purpose of this procedure is to simplify the analysis and facilitate comparison of results.
Conclusion
This study aimed to explore a novel method for obtaining long-term phreatic groundwater levels, by combining a data-driven time series model using transfer function-noise modelling with detrended historical meteorological time series that represent the current climate. In applying the method to an area in the Netherlands, we showed that groundwater observations can provide only a limited characterisation of drought events, especially in terms of extreme groundwater levels, due to their limited length. We conclude that additional insight into low groundwater levels can be provided by using a long-term detrended meteorological series and time series modelling.
A profound practical implication of our study is that proposed interventions to mitigate (the consequences of) groundwater droughts may differ significantly depending on the characterisation method chosen. Insights from our study can therefore help water managers to consider more comprehensive return period-based criteria for interventions such as irrigation bans or the placement of more structural measures including water system adaptations.
Appendix
The transfer function-noise model parameters that were used in this study are shown in Table A1.
Table A1. Model parameters
Open access
