{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4. Confined Aquifer Test - Nevada"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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": {},
"source": [
"This example is taken from Kruseman and de Ridder (1990). It is based on a pumping test conducted in a fractured tertiary volcanic aquifer near Yucca Mountains, Nevada, US. \n",
"\n",
"The borehole extends to 1219 m depth and penetrates a 400 m thick sequence of fractured volcanic rocks with water entry points. At every entry point, the head is more or less the same, so they have good hydraulic connection. The water table is about 470 m below land surface, which indicates confined conditions. Drawdown data was monitored at the well, named ***UE-25b#1*** and at an observation well, named ***UE-25a#1***, 110 m away. Three pumping tests were conducted at the site and will be the object of our investigation. \n",
"\n",
"Although conventional Darcy flow does not apply to flow in discrete fractures, the fracture network at this site is sufficiently pervasive to justify the application of Darcy’s law at the aquifer scale. Groundwater flow to the well occurs through fractures, while the matrix exchanges water with the fracture network. For implementation in `timflow`, the fractured aquifer system is approximated using a ModelMaq configuration. This is achieved by representing the matrix as an upper aquifer layer, separated by a 1 m thick aquitard from a lower aquifer layer representing the fracture network."
]
},
{
"cell_type": "markdown",
"metadata": {
"scrolled": true,
"tags": [
"hide-input"
]
},
"source": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" optimal \n",
" \n",
" \n",
" \n",
" \n",
" kaq_1_1 \n",
" 0.877139 \n",
" \n",
" \n",
" Saq_1_1 \n",
" 0.000005 \n",
" \n",
" \n",
" c_1_1 \n",
" 12.912907 \n",
" \n",
" \n",
" rc \n",
" 0.105678 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"hm1 = w.headinside(t1)\n",
"hm2 = ml.head(r, 0, t2)\n",
"plt.figure(figsize=(8, 5))\n",
"plt.semilogx(t1, -h1, \".\", label=\"obs UE-25b#1\")\n",
"plt.semilogx(t1, -hm1[-1], label=\"timflow UE-25b#1\")\n",
"plt.semilogx(t2, -h2, \".\", label=\"obs UE-25a#1\")\n",
"plt.semilogx(t2, -hm2[-1], label=\"timflow UE-25a#1\")\n",
"plt.xlabel(\"time [d]\")\n",
"plt.ylabel(\"drawdown [m]\")\n",
"plt.title(\"Model Results\")\n",
"plt.legend()\n",
"plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Comparison of results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The performance of `timflow` was evaluated by comparison to the results obtained from the software AQTESOLV (Duffield, 2007), MLU (Hemker & Post, 2014), and Kruseman and de Ridder (1990), here abbreviated to K&dR.\n",
"\n",
"The problem is solved in AQTESOLV using a double-porosity analytical solution proposed by Moench (1984). The MLU (Hemker & Post, 2014) solution used a similar approach to our `timflow` model by simulating the matrix as a very-low transmissivity aquifer on top of the fractured aquifer and separated by a zero-storage aquitard. However, the calibrated MLU model only fits the observations of well UE-25b#1. When the observations of well UE-25a#1 are included, the fit is not as good as initially thought. Kruseman et al. (1990) solved the problem using a graphical method, where the transmissivity was calculated as one aquifer and the storativity was separated between matrix and fractures. \n",
"\n",
"Overall, `timflow` showed similar results to MLU but with a better fit since both observation wells are taken into account with the `timflow` model. While the AQTESOLV obtained the best fit using Moench's analytical solution."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" \n",
" km [m/d] \n",
" Sm [1/m] \n",
" kf [m/d] \n",
" Sf [1/m] \n",
" c [d] \n",
" rc [m] \n",
" RMSE [m] \n",
" \n",
" \n",
" \n",
" \n",
" timflow \n",
" 0.00 \n",
" 3.75e-04 \n",
" 0.88 \n",
" 5.07e-06 \n",
" 12.91 \n",
" 0.11 \n",
" 0.198 \n",
" \n",
" \n",
" AQTESOLV \n",
" 0.15 \n",
" 5.48e-04 \n",
" 0.94 \n",
" 1.95e-03 \n",
" - \n",
" 0.11 \n",
" - \n",
" \n",
" \n",
" MLU \n",
" 0.00 \n",
" 3.75e-04 \n",
" 0.84 \n",
" 8.05e-06 \n",
" 5.20 \n",
" - \n",
" - \n",
" \n",
" \n",
" K&dR \n",
" 0.83 \n",
" 3.75e-04 \n",
" 0.83 \n",
" 4.00e-06 \n",
" - \n",
" - \n",
" - \n",
" \n",
" \n",
"
\n"
],
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = pd.DataFrame(\n",
" columns=[\n",
" \"km [m/d]\",\n",
" \"Sm [1/m]\",\n",
" \"kf [m/d]\",\n",
" \"Sf [1/m]\",\n",
" \"c [d]\",\n",
" \"rc [m]\",\n",
" \"RMSE [m]\",\n",
" ],\n",
" index=[\"timflow\", \"AQTESOLV\", \"MLU\", \"K&dR\"],\n",
")\n",
"\n",
"t.loc[\"timflow\"] = np.concatenate(\n",
" [np.array([0.00025, 3.750e-04]), cal.parameters[\"optimal\"].values, [cal.rmse()]]\n",
")\n",
"t.loc[\"AQTESOLV\"] = [0.1513, 5.4800e-4, 0.9366, 1.9508e-3, \"-\", 0.11, \"-\"]\n",
"t.loc[\"MLU\"] = [0.00025, 3.750e-04, 0.8351125, 8.05e-6, 5.2035, \"-\", \"-\"]\n",
"t.loc[\"K&dR\"] = [0.8325, 3.750e-4, 0.8325, 4.000e-6, \"-\", \"-\", \"-\"]\n",
"\n",
"t_formatted = t.style.format(\n",
" {\n",
" \"km [m/d]\": \"{:.2f}\",\n",
" \"Sm [1/m]\": \"{:.2e}\",\n",
" \"kf [m/d]\": \"{:.2f}\",\n",
" \"Sf [1/m]\": \"{:.2e}\",\n",
" \"c [d]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.2f}\",\n",
" \"rc [m]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.2f}\",\n",
" \"RMSE [m]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.3f}\",\n",
" }\n",
")\n",
"t_formatted"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reference\n",
"\n",
"* Duffield, G.M. (2007), AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n",
"* Hemker, K. en Post V. (2014) MLU for Windows: well flow modeling in multilayer aquifer systems; MLU User's guide. https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmicrofem.com%2Fdownload%2Fmlu-user.pdf&data=05%7C02%7CMark.Bakker%40tudelft.nl%7Cad7f16364d2d4fd55dbf08de73832eaa%7C096e524d692940308cd38ab42de0887b%7C0%7C0%7C639075204580287861%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=OBoe8seXZUfoat89Dfr4g6lF%2Bn1FdtXqtp%2F18BMXCn0%3D&reserved=0\n",
"* Kruseman, G.P. and De Ridder, N.A. (1990), Analysis and evaluation of pumping test data, 2nd edn. ILRI Publ. 47, ILRI, Wageningen, The Netherlands\n",
"* Moench, A.F. (1984), Double-Porosity Models for a Fissured Groundwater Reservoir With Fracture Skin, Water Resour. Res., 20( 7), 831– 846, doi:10.1029/WR020i007p00831."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "timflow (3.13.11)",
"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.11"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
![](\"figs/Nevada.png\")
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load Data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# time and drawdown of pumped well UE-25b#1\n",
"data1 = np.loadtxt(\"data/double-porosity-pumpingwell.txt\", skiprows=1)\n",
"t1 = data1[:, 0] # days\n",
"h1 = data1[:, 1] # m\n",
"\n",
"# time and drawdown of observation well UE-25a#1\n",
"data2 = np.loadtxt(\"data/double-porosity-110m.txt\", skiprows=1)\n",
"t2 = data2[:, 0]\n",
"h2 = data2[:, 1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Paramaters and model"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# known parameters\n",
"H = 400 # aquifer thickness in m\n",
"Q = 3093.12 # constant pumping rate in m^3/d\n",
"r = 110 # distance from pumping well to observation well in m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the `timflow` model, we will adopt the parameters for the first layer from the results obtained from Kruseman and de Ridder (1970). "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# parameters calculated by Kruseman and de Ridder (1970)\n",
"km = 0.1 / H # hydraulic conductivity of matrix\n",
"Sm = 3.75e-4 # specific storage of matrix"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"self.neq 1\n",
"solution complete\n"
]
}
],
"source": [
"ml = tft.ModelMaq(\n",
" kaq=[km, 1],\n",
" z=[0, -400, -401, -801],\n",
" c=5,\n",
" Saq=[Sm, 1e-3],\n",
" Sll=0,\n",
" topboundary=\"conf\",\n",
" tmin=1e-5,\n",
" tmax=3,\n",
")\n",
"w = tft.Well(ml, xw=0, yw=0, rw=0.11, rc=0, tsandQ=[0, Q], layers=1)\n",
"ml.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Estimate aquifer parameters\n",
"The aquifer parameters are estimated using head observations at both wells."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"..........................................................................................................................................................................................................\n",
"Fit succeeded.\n"
]
}
],
"source": [
"cal = tft.Calibrate(ml)\n",
"cal.set_parameter(name=\"kaq\", initial=10, layers=1)\n",
"cal.set_parameter(name=\"Saq\", initial=1e-4, pmin=0, layers=1)\n",
"cal.set_parameter(name=\"c\", initial=10, layers=1)\n",
"cal.set_parameter_by_reference(name=\"rc\", parameter=w.rc, initial=0)\n",
"\n",
"cal.seriesinwell(name=\"UE-25b#1\", element=w, t=t1, h=h1)\n",
"cal.series(name=\"UE-25a#1\", x=r, y=0, t=t2, h=h2, layer=1)\n",
"cal.fit(report=False)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"