{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1.Test for an anisotropic water-table aquifer - Moench Example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Import packages"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"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": {
"tags": []
},
"source": [
"### Introduction and Conceptual Model\n",
"\n",
"This test is based on a synthetic example reported by Barlow and Moench (1999), utilizing an analytical solution developed by Moench and Allen (1997) for the transient flow of partially-penetrating wells in unconfined aquifers. \n",
"\n",
"The aquifer is partially saturated with water (10 m water table). A pumping well is screened from 5 to 10 m depth. The well and the well-casing radius is 0.1 m. \n",
"\n",
"Drawdown is recorded at the pumping well and four piezometers located at two different distances and two different depths. Two piezometers, PS1 and PS2, are located at one-meter depth below the water table and 3.16 and 31.6 m distance, respectively. Another two (PD1 and PD2) piezometers are at 7.5 m depth below the water table and the same distances, directly below the previous piezometers. The figure below shows the location of the well and the piezometers\n",
"\n",
"Pumping starts at time t = 0 at a constant rate of 172.8 m3/d. Drawdown is recorded until t = 3 days."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"![](\"./figs/Moench.png\")
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load data\n",
"\n",
"The dataset for each well consists of a column with the time data in seconds and drawdown in meters. We are loading it and converting it to days and meters."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"data0 = np.loadtxt(\"data/moench_pumped.txt\", skiprows=1)\n",
"t0 = data0[:, 0] / 60 / 60 / 24 # convert time from seconds to days\n",
"h0 = -data0[:, 1] # converting drawdown to heads\n",
"data1 = np.loadtxt(\"data/moench_ps1.txt\", skiprows=1)\n",
"t1 = data1[:, 0] / 60 / 60 / 24 # convert time from seconds to days\n",
"h1 = -data1[:, 1]\n",
"data2 = np.loadtxt(\"data/moench_pd1.txt\", skiprows=1)\n",
"t2 = data2[:, 0] / 60 / 60 / 24 # convert time from seconds to days\n",
"h2 = -data2[:, 1]\n",
"data3 = np.loadtxt(\"data/moench_ps2.txt\", skiprows=1)\n",
"t3 = data3[:, 0] / 60 / 60 / 24 # convert time from seconds to days\n",
"h3 = -data3[:, 1]\n",
"data4 = np.loadtxt(\"data/moench_pd2.txt\", skiprows=1)\n",
"t4 = data4[:, 0] / 60 / 60 / 24 # convert time from seconds to days\n",
"h4 = -data4[:, 1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parameters and model"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"b = 10 # aquifer thickness, m\n",
"Q = 172.8 # constant discharge rate, m^3/d\n",
"rw = 0.1 # well radius, m\n",
"rc = 0.1 # casing radius, m\n",
"r1 = 3.16 # distance of closer wells, m\n",
"r2 = 31.6 # distance of wells more far away, m"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"self.neq 1\n",
"solution complete\n"
]
}
],
"source": [
"ml = tft.Model3D(\n",
" kaq=1,\n",
" z=[0, -0.1, -2.1, -5.1, -10.1],\n",
" Saq=[0.1, 1e-4, 1e-4, 1e-4],\n",
" kzoverkh=1,\n",
" tmin=1e-5,\n",
" tmax=3,\n",
" phreatictop=True,\n",
")\n",
"w = tft.Well(ml, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, Q)], layers=3)\n",
"ml.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Estimate aquifer parameters"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"................................................................\n",
"Fit succeeded.\n"
]
}
],
"source": [
"cal = tft.Calibrate(ml)\n",
"cal.set_parameter(name=\"kaq\", initial=1, layers=[0, 1, 2, 3])\n",
"cal.set_parameter(name=\"Saq\", initial=0.2, layers=[0])\n",
"cal.set_parameter(name=\"Saq\", initial=1e-4, pmin=0, layers=[1, 2, 3])\n",
"cal.set_parameter(name=\"kzoverkh\", initial=0.1, pmin=0, layers=[0, 1, 2, 3])\n",
"\n",
"cal.seriesinwell(name=\"pumped\", element=w, t=t0, h=h0)\n",
"\n",
"cal.series(name=\"PS1\", x=r1, y=0, t=t1, h=h1, layer=1)\n",
"cal.series(name=\"PD1\", x=r1, y=0, t=t2, h=h2, layer=3)\n",
"cal.series(name=\"PS2\", x=r2, y=0, t=t3, h=h3, layer=1)\n",
"cal.series(name=\"PD2\", x=r2, y=0, t=t4, h=h4, layer=3)\n",
"\n",
"cal.fit(report=False)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" optimal | \n",
"
\n",
" \n",
" \n",
" \n",
" | kaq_0_3 | \n",
" 9.067649 | \n",
"
\n",
" \n",
" | Saq_0_0 | \n",
" 0.172882 | \n",
"
\n",
" \n",
" | Saq_1_3 | \n",
" 0.000039 | \n",
"
\n",
" \n",
" | kzoverkh_0_3 | \n",
" 0.535051 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"hm0 = ml.head(0, 0, t0, layers=3)[0]\n",
"hm1 = ml.head(r1, 0, t1, layers=1)[0]\n",
"hm2 = ml.head(r1, 0, t2, layers=3)[0]\n",
"hm3 = ml.head(r2, 0, t3, layers=1)[0]\n",
"hm4 = ml.head(r2, 0, t4, layers=3)[0]\n",
"\n",
"plt.semilogx(t0, -h0, \".\", label=\"pumped\")\n",
"plt.semilogx(t0, -hm0, label=\"ttim pumped\")\n",
"plt.semilogx(t1, -h1, \".\", label=\"PS1\")\n",
"plt.semilogx(t1, -hm1, label=\"ttim PS1\")\n",
"plt.semilogx(t2, -h2, \".\", label=\"PD1\")\n",
"plt.semilogx(t2, -hm2, label=\"ttim PD1\")\n",
"plt.semilogx(t3, -h3, \".\", label=\"PS2\")\n",
"plt.semilogx(t3, -hm3, label=\"ttim PS2\")\n",
"plt.semilogx(t4, -h4, \".\", label=\"PD2\")\n",
"plt.semilogx(t4, -hm4, label=\"ttim PD2\")\n",
"\n",
"plt.xlabel(\"time (d)\")\n",
"plt.ylabel(\"drawdown (m)\")\n",
"plt.title(\"Model Results\")\n",
"plt.legend(bbox_to_anchor=(1.05, 1))\n",
"plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Comparison of results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The performance of `timflow` is compared with the stimulated results using Moench's parameters (Barlow and Moench, 1999). The RMSE decreased with the performance of `timflow`. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" | | \n",
" k [m/d] | \n",
" Sy [-] | \n",
" Ss [1/m] | \n",
" kz/kh | \n",
" RMSE [m] | \n",
"
\n",
" \n",
" \n",
" \n",
" | timflow | \n",
" 9.07 | \n",
" 0.2 | \n",
" 3.87e-05 | \n",
" 0.5 | \n",
" 0.0104 | \n",
"
\n",
" \n",
" | Moench | \n",
" 8.64 | \n",
" 0.2 | \n",
" 2.00e-05 | \n",
" 0.5 | \n",
" 0.0613 | \n",
"
\n",
" \n",
"
\n"
],
"text/plain": [
""
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = pd.DataFrame(\n",
" columns=[\"k [m/d]\", \"Sy [-]\", \"Ss [1/m]\", \"kz/kh\", \"RMSE [m]\"],\n",
" index=[\"timflow\", \"Moench\"],\n",
")\n",
"\n",
"t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
"t.loc[\"Moench\"] = [8.64, 0.2, 2e-5, 0.5, 0.061318]\n",
"\n",
"t_formatted = t.style.format(\n",
" {\n",
" \"k [m/d]\": \"{:.2f}\",\n",
" \"Sy [-]\": \"{:.1f}\",\n",
" \"Ss [1/m]\": \"{:.2e}\",\n",
" \"kz/kh\": \"{:.1f}\",\n",
" \"RMSE [m]\": \"{:.4f}\",\n",
" }\n",
")\n",
"t_formatted"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"* Barlow, P.M., Moench, A.F., 1999. WTAQ, a computer program for calculating drawdowns and estimating hydraulic properties for confined and water-table aquifers. 99-4225, US Dept. of the Interior, US Geological Survey\n",
"* Moench, Allen, F., 1997. Flow to a well of finite diameter in a homogeneous, anisotropic water table aquifer. Water Resources Research 34, 2431–2432."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.13.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}