{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Standardized Groundwater Index\n",
"*R.A. Collenteur, University of Graz, November 2020*\n",
"\n",
"To study the occurrence of groundwater droughts [Bloomfield and Marchant (2013)](#References) developed the Standardized Groundwater Index (SGI). More and more SGI values are being used to study and quantify groundwater droughts. In this Notebook it is shown how to compute the SGI using Pastas, and how Pastas models may be used to obtain groundwater level time series with regular time steps. The SGI implemented in Pastas (`ps.stats.sgi`) is based on the description in [Bloomfield and Marchant (2013)](#References).\n",
"\n",
"The SGI requires regular time steps between groundwater levels observation, while historic groundwater level time series are often characterized by irregular time intervals between observations. To overcome this issue, [Marchant and Bloomfield(2018)](#References) applied time series models using impulse response functions to simulate groundwater level time series at a regular time interval. Here, this methodology is extended by using evaporation and precipitation as model input and using a nonlinear recharge model ([Collenteur et al. (2021)](#References)) to compute groundwater recharge and finally groundwater levels.\n",
"\n",
"**Note that this notebook is meant as an example of how Pastas models may be used to support studies computing SGI values, and not as an guide how to compute or interpret the SGI values.**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import pastas as ps\n",
"import matplotlib.pyplot as plt\n",
"\n",
"ps.set_log_level(\"ERROR\")\n",
"ps.show_versions()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The first example model\n",
"\n",
"### 1. Loading the data\n",
"In this example we model the groundwater levels for a monitoring well (B32C0639, filter 1) near the town \"de Bilt\" in the Netherlands. Precipitation and evaporation are available from the nearby meteorological station of the KNMI. The groundwater level observations have irregular time steps."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load input data\n",
"head = pd.read_csv(\n",
" \"data/B32C0639001.csv\", parse_dates=[\"date\"], index_col=\"date\"\n",
").squeeze()\n",
"evap = pd.read_csv(\"data/evap_260.csv\", index_col=0, parse_dates=[0]).squeeze()\n",
"rain = pd.read_csv(\"data/rain_260.csv\", index_col=0, parse_dates=[0]).squeeze()\n",
"\n",
"# Plot input data\n",
"ps.plots.series(head, [rain, evap]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. Creating and calibrating the model\n",
"We now create a simple time series model using a parsimonious non-linear recharge model to translate precipitation and evaporation into groundwater recharge. The recharge flux is then convolved with an exponential response function to compute the contribution of the recharge to the groundwater level fluctuations. The results are plotted below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create the basic Pastas model\n",
"ml = ps.Model(head)\n",
"ml.add_noisemodel(ps.ArNoiseModel())\n",
"\n",
"# Add a recharge model\n",
"rch = ps.rch.FlexModel()\n",
"rm = ps.RechargeModel(rain, evap, recharge=rch, rfunc=ps.Exponential(), name=\"rch\")\n",
"ml.add_stressmodel(rm)\n",
"\n",
"# Solve the model\n",
"ml.solve(tmin=\"1990\", report=False)\n",
"ml.plots.results(figsize=(10, 6));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. Computing and visualizing the SGI\n",
"The plot above shows that we have a pretty good model fit with the data. This is particularly important when we want to compute the SGI using the simulated time series. We now compute the SGI and show the models results and estimated SGI in one figure. A possible extension to the SGI computation below is to take the uncertainty of the groundwater level simulation into account, as is done by [Marchant and Bloomfield (2018)](#References)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute the SGI\n",
"sim = ml.simulate(tmin=\"1990\")\n",
"sgi = ps.stats.sgi(sim.resample(\"W\").mean())\n",
"ci = ml.solver.prediction_interval(n=10)\n",
"\n",
"# Make the plot\n",
"fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(10, 5), sharex=True)\n",
"\n",
"# Upper subplot\n",
"sim.plot(ax=ax1, zorder=10)\n",
"ml.oseries.series.plot(ax=ax1, linestyle=\" \", marker=\".\", color=\"k\")\n",
"ax1.fill_between(ci.index, ci.iloc[:, 0], ci.iloc[:, 1], color=\"gray\")\n",
"ax1.legend([\"Simulation\", \"Observations\", \"Prediction interval\"], ncol=3)\n",
"\n",
"# Lower subplot\n",
"sgi.plot(ax=ax2, color=\"k\")\n",
"ax2.axhline(0, linestyle=\"--\", color=\"k\")\n",
"droughts = sgi.to_numpy(copy=True)\n",
"droughts[droughts > 0] = 0\n",
"ax2.fill_between(sgi.index, 0, droughts, color=\"C0\")\n",
"\n",
"# Dress up the plot\n",
"ax1.set_ylabel(\"GWL [m]\")\n",
"ax1.set_title(\"Groundwater levels\")\n",
"ax2.set_ylabel(\"SGI [-]\")\n",
"ax2.set_title(\"Standardized Groundwater Index\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Second Example: data with trends\n",
"For the example above precipitation and evaporation were sufficient to accurately simulate the groundwater levels. Now we look at an example of where this is not the case. The groundwater levels are again observed near the town of de Bilt in the Netherlands. The time series have a more irregularities in the time step between observations and end with high frequency observations.\n",
"\n",
"### 1. Create a simple model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Loads heads and create Pastas model\n",
"head2 = pd.read_csv(\"data/B32C0609001.csv\", parse_dates=[0], index_col=0).squeeze()\n",
"ml2 = ps.Model(head2)\n",
"ml2.add_noisemodel(ps.ArNoiseModel())\n",
"\n",
"# Add a recharge model\n",
"rch = ps.rch.FlexModel()\n",
"rm = ps.RechargeModel(rain, evap, recharge=rch, rfunc=ps.Exponential(), name=\"rch\")\n",
"ml2.add_stressmodel(rm)\n",
"\n",
"# Solve and plot the model\n",
"ml2.solve(tmin=\"1990\", report=False)\n",
"ml2.plots.results(figsize=(10, 6));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. Add linear trend\n",
"Clearly the model fit with the data in the above figure is not so good. Looking at the model residuals (simulation - observation) we can observe a steady upward trend in the residuals. Let's try and add a linear trend to the model to improve the groundwater level simulation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Add a linear trend\n",
"tm = ps.LinearTrend(\"1990\", \"2020\", name=\"trend\")\n",
"ml2.add_stressmodel(tm)\n",
"\n",
"# Solve the model\n",
"ml2.del_noisemodel()\n",
"# ml2.solve(tmin=\"1990\", report=False) # Get better initial estimated first\n",
"ml2.add_noisemodel(ps.ArNoiseModel())\n",
"ml2.solve(tmin=\"1990\", report=False)\n",
"ml2.plots.results(figsize=(10, 6));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. Computing and plotting the SGI\n",
"The model fit for the model above looks a lot better. Now we can compute and plot the SGI again as we did before."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute the SGI\n",
"sim = ml2.simulate(tmin=\"1990\")\n",
"sgi = ps.stats.sgi(sim.resample(\"W\").mean())\n",
"ci = ml2.solver.prediction_interval(n=10)\n",
"\n",
"# Make the plot\n",
"fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(10, 5), sharex=True)\n",
"\n",
"# Upper subplot\n",
"sim.plot(ax=ax1, zorder=10)\n",
"ml2.oseries.series.plot(ax=ax1, linestyle=\" \", marker=\".\", color=\"k\")\n",
"ax1.fill_between(ci.index, ci.iloc[:, 0], ci.iloc[:, 1], color=\"gray\")\n",
"ax1.legend([\"Simulation\", \"Observations\", \"Prediction interval\"], ncol=3)\n",
"\n",
"# Lower subplot\n",
"sgi.plot(ax=ax2, color=\"k\")\n",
"ax2.axhline(0, linestyle=\"--\", color=\"k\")\n",
"droughts = sgi.to_numpy(copy=True)\n",
"droughts[droughts > 0] = 0\n",
"ax2.fill_between(sgi.index, 0, droughts, color=\"C0\")\n",
"\n",
"# Dress up the plot\n",
"ax1.set_ylabel(\"GWL [m]\")\n",
"ax1.set_title(\"Groundwater levels\")\n",
"ax2.set_ylabel(\"SGI [-]\")\n",
"ax2.set_title(\"Standardized Groundwater Index\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What about human influenced groundwater systems?\n",
"Let's explore the possibilities of using the Pastas framework a bit here. The first example showed SGI values for a system under natural conditions, with only recharge being enough to explain the groundwater level fluctuations. In the second example a small linear trend had to be added, without explicit knowledge of what may have caused this trend. In this third and final example we consider an aquifer system that is influenced by groundwater pumping. \n",
"\n",
"The question we want to answer is how the SGI values may have looked without groundwater pumping (a natural system) and compare these to the SGI values with groundwater pumping. We can see clearly from the model that groundwater pumping decreased the groundwater levels, but how does it impact the SGI values?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load input data\n",
"head = pd.read_csv(\"data_notebook_9/head.csv\", parse_dates=True, index_col=0).squeeze()\n",
"prec = pd.read_csv(\"data_notebook_9/prec.csv\", parse_dates=True, index_col=0).squeeze()\n",
"evap = pd.read_csv(\"data_notebook_9/evap.csv\", parse_dates=True, index_col=0).squeeze()\n",
"well = pd.read_csv(\"data_notebook_9/well.csv\", parse_dates=True, index_col=0).squeeze()\n",
"\n",
"# ps.validate_stress(well)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the well data is not equidistant. For this example, we apply a backfill after resampling the time series to daily values."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"well = well.asfreq(\"D\").bfill()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are ready to build the model."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create the Pastas model\n",
"ml3 = ps.Model(head, name=\"heads\")\n",
"ml3.add_noisemodel(ps.ArNoiseModel())\n",
"\n",
"# Add recharge and a well\n",
"sm = ps.RechargeModel(\n",
" prec, evap, ps.Exponential(), name=\"rch\", recharge=ps.rch.FlexModel()\n",
")\n",
"wm = ps.StressModel(well, ps.Exponential(), well.name, up=False, settings=\"well\")\n",
"ml3.add_stressmodel([sm, wm])\n",
"\n",
"# Solve the model\n",
"ml3.solve(report=False)\n",
"ml3.plots.results(figsize=(10, 6));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. SGI with and without groundwater pumping\n",
"Now that we have a model with a reasonably good fit, we can use the model to separate the effect of groundwater pumping from the effect of recharge. We then compute SGI values on the groundwater levels with and without pumping and compare them visually. The results are shown below, and show very different SGI values as expected."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute the SGI\n",
"sim = ml3.simulate(tmin=\"1940\")\n",
"sgi = ps.stats.sgi(sim.resample(\"M\").mean())\n",
"recharge = ml3.get_contribution(\"rch\", tmin=\"1940\")\n",
"sgi2 = ps.stats.sgi(recharge.resample(\"M\").mean())\n",
"# ci = ml3.solver.prediction_interval()\n",
"\n",
"# Make the plot\n",
"fig, [ax1, ax2, ax3] = plt.subplots(3, 1, figsize=(10, 6), sharex=True)\n",
"\n",
"sim.plot(ax=ax1, x_compat=True)\n",
"(recharge + ml3.get_parameters(\"constant\")).plot(ax=ax1, linestyle=\"--\")\n",
"ml3.oseries.series.plot(\n",
" ax=ax1, linestyle=\" \", marker=\".\", zorder=-1, markersize=2, color=\"k\", x_compat=True\n",
")\n",
"# ax1.fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], color=\"gray\")\n",
"ax1.legend([\"Simulation\", \"Simulation w/o pumping\"], ncol=1)\n",
"\n",
"sgi.plot(ax=ax2, color=\"k\", x_compat=True)\n",
"ax2.axhline(0, linestyle=\"--\", color=\"k\")\n",
"droughts = sgi.to_numpy(copy=True)\n",
"droughts[droughts > 0] = 0\n",
"ax2.fill_between(sgi.index, 0, droughts, color=\"C0\")\n",
"\n",
"sgi2.plot(ax=ax3, color=\"k\", x_compat=True)\n",
"ax3.axhline(0, linestyle=\"--\", color=\"k\")\n",
"droughts = sgi2.to_numpy(copy=True)\n",
"droughts[droughts > 0] = 0\n",
"ax3.fill_between(sgi2.index, 0, droughts, color=\"C1\")\n",
"\n",
"ax1.set_ylabel(\"GWL [m]\")\n",
"ax1.set_title(\"Groundwater levels\")\n",
"ax2.set_ylabel(\"SGI [-]\")\n",
"ax2.set_title(\"SGI With Groundwater pumping\")\n",
"ax3.set_ylabel(\"SGI [-]\")\n",
"ax3.set_title(\"SGI under 'Natural conditions'\")\n",
"plt.xlim(pd.Timestamp(\"1940\"), pd.Timestamp(\"2016\"));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"- Bloomfield, J. P. and Marchant, B. P.: [Analysis of groundwater drought building on the standardised precipitation index approach](https://hess.copernicus.org/articles/17/4769/2013/), Hydrol. Earth Syst. Sci., 17, 4769–4787, 2013.\n",
"- Marchant, B. and Bloomfield, J.: [Spatio-temporal modelling of the status of groundwater droughts](https://doi.org/10.1016/j.jhydrol.2018.07.009), J. Hydrol., 564, 397–413, 2018\n",
"- Collenteur, R., 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](https://doi.org/10.5194/hess-2020-392), Hydrol. Earth Syst. Sci., 25, 2931–2949.\n",
"\n",
"## Data Sources\n",
"\n",
"- The precipitation and evaporation time series are taken from the Dutch KNMI, meteorological station \"de Bilt\" (www.knmi.nl).\n",
"- The groundwater level time series were downloaded from Dinoloket (www.dinoloket.nl)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
},
"vscode": {
"interpreter": {
"hash": "29475f8be425919747d373d827cb41e481e140756dd3c75aa328bf3399a0138e"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}