{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# 2. Unconfined Aquifer Test - Vennebulten"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Import packages"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import timflow.transient as tft\n",
"\n",
"plt.rcParams[\"figure.figsize\"] = (5, 3) # default figure size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Introduction and Conceptual Model"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"In aquifer tests in unconfined aquifers, there is also the vertical component to flow to the well. The drawdown data shows the delayed water table response, a distinguishable S-shape in the log-log plot. In the early times of the drawdown, the drawdown behaves as a confined aquifer: when the aquifer releases the elastic storage. However, as pumping continues, the water table storage begins to be released, generating further drawdown and the S-shape.\n",
"\n",
"This test conducted in Vennebulten, the Netherlands, is reported in Kruseman et al. (1970). The cross-section consists of a first layer up to 6 m depth of very fine and loamy sands, followed by coarse sands until 21 m deep.\n",
"\n",
"The screen of the pumping well is placed between 10 and 21 meters depth, and pumping has taken place for 25 hours at a rate of 873 m3/d. The available drawdown data comes from two piezometers, a shallow one, screened at 3 m depth, and a deeper one, screened in the depths between 12 to 19 m. Both wells are located 90 m from the pumping well."
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" optimal \n",
" \n",
" \n",
" \n",
" \n",
" Saq_1_11 \n",
" 0.000140 \n",
" \n",
" \n",
" kaq_0_11 \n",
" 0.080617 \n",
" \n",
" \n",
" Saq_12_26 \n",
" 0.000029 \n",
" \n",
" \n",
" kaq_12_26 \n",
" 121.610702 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"hs_1 = ml.head(r, 0, ts, layers=3)\n",
"hd_1 = wobs.headinside(td)\n",
"\n",
"plt.semilogx(ts, hs, \".\", label=\"shallow obs\")\n",
"plt.semilogx(ts, hs_1[0], label=\"shallow ttim\")\n",
"plt.semilogx(td, hd, \".\", label=\"deep obs\")\n",
"plt.semilogx(td, hd_1[0], \"--\", label=\"deep ttim\")\n",
"plt.xlabel(\"time [d]\")\n",
"plt.ylabel(\"head change [m]\")\n",
"plt.title(\"Timflow Unconfined Model Results - Shallow Piezometer\")\n",
"plt.legend()\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Comparison of results\n",
"The performance of `timflow` is compared with a model performed with MLU (Hemker & Post, 2014). The schematization of the MLU model differs considerably. The model consists of three layers: a very thin fine-grained upper aquifer layer with horizontal flow and two deeper sublayers with a thickness of 1 and 10 meters. In total, four parameters were estimated (Ss0, c2, Ss2, k2), while fixed values were chosen for the remaining ones. \n",
"\n",
"The RMSE of MLU is smaller than that of timflow, but MLU uses more layers with fixed parameters."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" \n",
" Ss0 [1/m] \n",
" k0 [m/d] \n",
" c1 [d] \n",
" Ss1 [1/m] \n",
" k1 [m/d] \n",
" c2 [d] \n",
" Ss2 [1/m] \n",
" k2 [m/d] \n",
" RMSE [m] \n",
" \n",
" \n",
" \n",
" \n",
" timflow \n",
" 1.40e-04 \n",
" 0.08 \n",
" - \n",
" 2.94e-05 \n",
" 121.61 \n",
" - \n",
" - \n",
" - \n",
" 0.0082 \n",
" \n",
" \n",
" MLU \n",
" 3.80e-01 \n",
" 1.00 \n",
" 100 \n",
" 1.00e-06 \n",
" 50.00 \n",
" 209 \n",
" 5.58e-03 \n",
" 165.90 \n",
" 0.0005 \n",
" \n",
" \n",
"
\n"
],
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = pd.DataFrame(\n",
" columns=[\n",
" \"Ss0 [1/m]\",\n",
" \"k0 [m/d]\",\n",
" \"c1 [d]\",\n",
" \"Ss1 [1/m]\",\n",
" \"k1 [m/d]\",\n",
" \"c2 [d]\",\n",
" \"Ss2 [1/m]\",\n",
" \"k2 [m/d]\",\n",
" \"RMSE [m]\",\n",
" ],\n",
" index=[\"timflow\", \"MLU\"],\n",
")\n",
"\n",
"t.loc[\"timflow\"] = [\n",
" cal.parameters[\"optimal\"].values[0],\n",
" cal.parameters[\"optimal\"].values[1],\n",
" \"-\",\n",
" cal.parameters[\"optimal\"].values[2],\n",
" cal.parameters[\"optimal\"].values[3],\n",
" \"-\",\n",
" \"-\",\n",
" \"-\",\n",
" cal.rmse(),\n",
"]\n",
"t.loc[\"MLU\"] = [0.38, 1, 100, 1e-6, 50, 209, 0.00558, 165.9, 0.0005]\n",
"\n",
"t_formatted = t.style.format(\n",
" {\n",
" \"Ss0 [1/m]\": \"{:.2e}\",\n",
" \"k0 [m/d]\": \"{:.2f}\",\n",
" \"c1 [d]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.0f}\",\n",
" \"Ss1 [1/m]\": \"{:.2e}\",\n",
" \"k1 [m/d]\": \"{:.2f}\",\n",
" \"c2 [d]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.0f}\",\n",
" \"Ss2 [1/m]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.2e}\",\n",
" \"k2 [m/d]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.2f}\",\n",
" \"RMSE [m]\": \"{:.4f}\",\n",
" }\n",
")\n",
"t_formatted"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## References\n",
"* Hemker, K. en Post V. (2014) MLU for Windows: well flow modeling in multilayer aquifer systems; [MLU User's guide](https://microfem.com/download/mlu-user.pdf)\n",
"* Kruseman, G.P., De Ridder, N.A., Verweij, J.M., 1970. Analysis and evaluationof pumping test data. volume 11. International institute for land reclamation and improvement The Netherlands.\n",
"* Neuman, S.P., Witherspoon, P.A., 1969. Applicability of current theories of flow in leaky aquifers. Water Resources Research 5, 817–829."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:base] *",
"language": "python",
"name": "conda-base-py"
},
"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.13.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
![](\"./figs/Vennebulten.png\")
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load data "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"data1 = np.loadtxt(\"data/venne_shallow.txt\", skiprows=1)\n",
"ts = data1[:, 0] / 60 / 24 # convert min to days\n",
"hs = data1[:, 1]\n",
"\n",
"data2 = np.loadtxt(\"data/venne_deep.txt\", skiprows=1)\n",
"td = data2[:, 0] / 60 / 24 # convert min to days\n",
"hd = data2[:, 1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Parameters and model"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"b = -21 # aquifer thickness, m\n",
"r = 90 # distance from observation wells to pumping well, m\n",
"Q = 873 # constant discharge, m^3/d\n",
"k = 12 * [0.5] + 15 * [150]\n",
"z = np.hstack((np.arange(0, -6, -0.5), np.arange(-6, -21.1, -1)))\n",
"Ss = [0.2] + 26 * [1e-4]\n",
"kzoverkh = 0.1"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"self.neq 20\n",
"solution complete\n"
]
}
],
"source": [
"ml = tft.Model3D(\n",
" kaq=k,\n",
" z=z,\n",
" Saq=Ss,\n",
" kzoverkh=kzoverkh,\n",
" topboundary=\"phreatic\",\n",
" tmin=1e-4,\n",
" tmax=1.1,\n",
")\n",
"w = tft.Well(ml, xw=0, yw=0, rw=0.1, tsandQ=[(0, Q)], layers=range(10, 21))\n",
"wobs = tft.Well(ml, xw=90, yw=0, rw=0.1, tsandQ=[(0, 0)], layers=range(18, 27))\n",
"ml.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Estimate aquifer parameters"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"...............................................................................\n",
"Fit succeeded.\n",
"[[Fit Statistics]]\n",
" # fitting method = leastsq\n",
" # function evals = 76\n",
" # data points = 48\n",
" # variables = 4\n",
" chi-square = 0.00319397\n",
" reduced chi-square = 7.2590e-05\n",
" Akaike info crit = -453.649150\n",
" Bayesian info crit = -446.164346\n",
"[[Variables]]\n",
" Saq_1_11: 1.3997e-04 (init = 0.0001)\n",
" kaq_0_11: 0.08061669 (init = 1)\n",
" Saq_12_26: 2.9428e-05 (init = 0.0001)\n",
" kaq_12_26: 121.610702 (init = 150)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"......."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"."
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
".\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fit succeeded.\n",
"[[Fit Statistics]]\n",
" # fitting method = leastsq\n",
" # function evals = 76\n",
" # data points = 48\n",
" # variables = 4\n",
" chi-square = 0.00319397\n",
" reduced chi-square = 7.2590e-05\n",
" Akaike info crit = -453.649150\n",
" Bayesian info crit = -446.164346\n",
"[[Variables]]\n",
" Saq_1_11: 1.3997e-04 (init = 0.0001)\n",
" kaq_0_11: 0.08061669 (init = 1)\n",
" Saq_12_26: 2.9428e-05 (init = 0.0001)\n",
" kaq_12_26: 121.610702 (init = 150)\n"
]
}
],
"source": [
"cal = tft.Calibrate(ml)\n",
"\n",
"cal.set_parameter(\n",
" name=\"Saq\", initial=1e-4, layers=list(range(1, 12)), pmin=1e-6, pmax=1e-2\n",
")\n",
"cal.set_parameter(name=\"kaq\", initial=150, layers=list(range(0, 12)), pmin=0.05, pmax=1)\n",
"\n",
"cal.set_parameter(\n",
" name=\"Saq\", initial=1e-4, layers=list(range(12, 27)), pmin=1e-6, pmax=1e-2\n",
")\n",
"cal.set_parameter(name=\"kaq\", initial=150, layers=list(range(12, 27)), pmin=100, pmax=200)\n",
"\n",
"cal.series(name=\"shallow_obs\", x=r, y=0, t=ts, h=hs, layer=5)\n",
"cal.seriesinwell(name=\"wobs\", element=wobs, t=td, h=hd)\n",
"\n",
"cal.fit(report=True)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"