{
"cells": [
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"hideCode": false,
"hidePrompt": false,
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from myst_nb import glue\n",
"%config InlineBackend.rc = {'figure.dpi': 300}\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import sys\n",
"sys.path.append(\"../py_figures\")\n",
"import diagrams_2d_chap9 as diagrams_2d\n",
"import inter_fig_chap9"
]
},
{
"cell_type": "markdown",
"metadata": {
"hideCode": false,
"hidePrompt": false
},
"source": [
"(content:chap9.4:wellbore_storage)=\n",
"# Wells with wellbore storage and skin effect\n",
"\n",
"The previous solutions for transient flow to a well were based on two approximations: The boundary condition at the well was approximated by taking the limit for the well radius going to zero (Eq. {eq}`eq_bctheis19.1`), and the water level in the well is approximated by the head at the well radius in the aquifer. The first approximation neglects the storage of the water inside the well, which is referred to as wellbore storage. The second approximation neglects the resistance of the well screen, which is referred to as the skin effect. Both are taken into account in this section.\n",
"
Consider radial flow to a well in a confined aquifer with transmissivity $T$ and storage coefficient $S$. The well starts pumping with discharge $Q$ at time $t = 0$. The radius of the well screen is $r_w$ (see {numref}`Figure %s `). The resistance to flow through the well screen is equal to $c_w$ so that the radial component of the discharge vector right at the well screen may be written as\n",
"\n",
"$$\n",
"Q_r|_{r=r_w} = H\\frac{h_w - h|_{r=r_w}}{c_w}\n",
"$$(eq_Q_r9.4)\n",
"\n",
"where $h_w$ is the water level inside the well. Conversely, the water level in the well may be written as\n",
"\n",
"$$\n",
"h_w = h|_{r=r_w} + \\frac{c_w}{H}Q_r|_{r=r_w} \n",
"$$(eq_wbs_eq19.4)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"hideCode": false,
"hideOutput": true,
"hidePrompt": false,
"tags": [
"remove-input",
"remove-output"
]
},
"outputs": [
{
"data": {
"application/papermill.record/image/svg+xml": "\n\n\n",
"application/papermill.record/text/plain": ""
},
"metadata": {
"scrapbook": {
"mime_prefix": "application/papermill.record/",
"name": "fig_wellbore"
}
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"wellbore = diagrams_2d.wellbore()\n",
"glue(\"fig_wellbore\", wellbore, display=False)"
]
},
{
"cell_type": "markdown",
"metadata": {
"hideCode": false,
"hideOutput": true,
"hidePrompt": false
},
"source": [
"```{glue:figure} fig_wellbore\n",
"---\n",
"name: fig_well_bore9.4\n",
"---\n",
"Vertical cross section through a well with wellbore storage and skin effect. The head in the aquifer at the well screen is higher than the water level $h_w$ in the well because of the skin effect.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"hideCode": false,
"hidePrompt": false
},
"source": [
" The water balance for water inside the well is\n",
"\n",
"$$\n",
"A_w \\frac{\\text{d}h_w}{\\text{d}t} = Q_\\text{in} - Q_\\text{out}\n",
"$$(eq_waterbalanceinwell9.4)\n",
"\n",
"where $A_w$ is the cross-sectional area of the well screen at the level of the water table in the well ($A_w=\\pi r_w^2$ for wells with a constant well radius). The inflow $Q_\\text{in}$ consists of the inflow integrated along the well screen\n",
"\n",
"$$\n",
"Q_\\text{in} = -2\\pi r_w Q_r|_{r=r_w}\n",
"$$(eq_wbs_eq29.4)\n",
"\n",
"while the outflow is equal to the discharge of the well \n",
"\n",
"$$\n",
"Q_\\text{out} = Q {\\hskip 2em} t\\ge 0\n",
"$$(eq_wbs_eq39.4)\n",
"\n",
"Substitution of Eq. ({eq}`eq_wbs_eq19.4`) for $h_w$, Eq. ({eq}`eq_wbs_eq29.4`) for $Q_\\text{in}$, and Eq. ({eq}`eq_wbs_eq39.4`) for $Q_\\text{out}$ in the water balance inside the well gives\n",
"\n",
"$$\n",
"A_w \\left(\\left.\\frac{\\partial h}{\\partial t}\\right|_{r=r_w} + \\frac{c_w}{H} \\left.\\frac{\\partial Q_r}{\\partial t}\\right|_{r=r_w}\\right) = -2\\pi r_w Q_{r}|_{r=r_w} - Q {\\hskip 2em} t\\ge 0\n",
"$$(eq_waterbalance_inside_well9.4)\n",
"\n",
"
The governing differential equation for flow in the aquifer is the diffusion equation ({eq}`eq_deq_rad_sym_transient_flow9.0`). The Laplace transform of the diffusion equation is the modified Helmholtz equation\n",
"\n",
"$$\n",
"\\frac{\\text{d}^2 \\overline{h}}{\\text{d} r^2}+\\frac{1}{r}\\frac{\\text{d}\\overline{h}}{\\text{d}r}= \\frac{Sp}{T}\\overline{h}=w\\overline{h}\n",
"$$(eq_diffeq_laptran9.4)\n",
"\n",
"where \n",
"\n",
"$$\n",
"w = Sp/T\n",
"$$(eq_w9.4)\n",
"\n",
"The Laplace transform of the water balance inside the well (Eq. {eq}`eq_waterbalance_inside_well9.4`) is\n",
"\n",
"$$\n",
"A_w p \\overline{h}|_{r=r_w} + \\left(\\frac{A_w c_w p}{H} + 2\\pi r_w\\right) \\overline{Q}_r|_{r=r_w} = - \\frac{Q}{p}\n",
"$$(eq_ltwatbal9.4)\n",
"\n",
"
The general solution to the Laplace-transformed diffusion equation in radial coordinates (Eq. {eq}`eq_diffeq_laptran9.4`), for the case that h equals zero far away, is\n",
"\n",
"$$\n",
"\\overline{h}=B\\text{K}_0(r\\sqrt{w})\n",
"$$(eq_h_over9.4)\n",
"\n",
"where $B$ is a parameter to be determined from the boundary conditions. The corresponding discharge vector in the Laplace domain is\n",
"\n",
"$$\n",
"\\overline{Q}_r=-T\\frac{\\text{d}\\overline{h}}{\\text{d}r} = BT\\sqrt{w}\\text{K}_1(r\\sqrt{w})\n",
"$$(eq_Q_r_over9.4)\n",
"\n",
"Substitution of the equations for $h$ and $Q_r$ in the Laplace-transformed water balance (Eq. {eq}`eq_ltwatbal9.4`) gives\n",
"\n",
"$$\n",
"B = -\\frac{Q}{p} \\frac{1}{A_w p \\text{K}_0(r_w\\sqrt{w}) + \\left(\\frac{A_w c_w p}{H} + 2\\pi r_w\\right) T \\sqrt{w} \\text{K}_1(r_w \\sqrt{w}) }\n",
"$$(eq_B_sub9.4)\n",
"\n",
"Back-transformation of the Laplace-transformed solution to the time domain is again carried out using the Stehfest algorithm (see {numref}`content:chap5:laplace_transforms`). Note that for the case that wellbore storage is neglected, the cross-sectional area $A_w$ is set to zero in the water balance equation ({eq}`eq_waterbalanceinwell9.4`), which means that $(\\pi r_w^2=0)$ and $B$ reduces to \n",
"\n",
"$$\n",
"B = -\\frac{Q}{p} \\frac{1}{2\\pi r_w T \\sqrt{w} \\text{K}_1(r_w \\sqrt{w})}\n",
"$$(eq_B2.94)\n",
"\n",
"which is the solution when the boundary condition at the well screen is\n",
"\n",
"$$\n",
"-2\\pi r_w Q_r|_{r=r_w} = Q {\\hskip 2em} t\\ge 0\n",
"$$(eq_bcwelltrue9.4)\n",
"\n",
"In the example below, the head change is compared for a well without wellbore storage (boundary condition {eq}`eq_bcwelltrue9.4`), a well with wellbore storage but no skin effect, and a well with both wellbore storage and skin effect ({numref}`Figure %s `). For the values of the example, wellbore storage and skin effect have an effect on the head during the first few hours of pumping ($t < 0.1$ d). For the case with skin effect, the head inside the well starts to differ from the head just outside the well after a few minutes of pumping ($t > 0.002$ d)."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"hideCode": false,
"hidePrompt": false
},
"outputs": [],
"source": [
"# parameters\n",
"T = 200 # transmissivity, m^2/d\n",
"S = 0.005 # storage coefficient of aquifer, -\n",
"H = 10 # aquifer thickness, m\n",
"Q = 800 # discharge of well, m^3/d\n",
"rw = 0.3 # well radius, m\n",
"cw = 0.02 # resistance of well screen, d"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"hideCode": false,
"hidePrompt": false,
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"from scipy.special import factorial as fac\n",
"def stehfest_coef(M):\n",
" assert M % 2 == 0, 'M must be even'\n",
" M2 = int(M / 2)\n",
" V = np.zeros(M)\n",
" for k in range(1, M + 1):\n",
" z = np.zeros(min(k, M2) + 1)\n",
" for j in range(int((k + 1) / 2), min(k, M2) + 1):\n",
" z[j] = j ** M2 * fac(2 * j) / \\\n",
" (fac(M2 - j) * fac(j) * fac(j - 1) * fac(k - j) * fac(2 * j - k))\n",
" V[k - 1] = (-1) ** (k + M2) * np.sum(z)\n",
" return V"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"hideCode": false,
"hidePrompt": false,
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"def stehfest(x, t, func, M=12, **kwargs):\n",
" t = np.atleast_1d(t)\n",
" x = np.atleast_1d(x)\n",
" f = np.zeros((len(x), len(t)))\n",
" V = stehfest_coef(M)\n",
" for i in range(len(t)):\n",
" p = np.arange(1, M + 1) * np.log(2) / t[i]\n",
" for j in range(len(x)):\n",
" fbar = func(p, x[j], **kwargs)\n",
" f[j, i] = np.log(2) / t[i] * np.sum(V * fbar)\n",
" return f"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"hideCode": false,
"hidePrompt": false,
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"from scipy.special import exp1\n",
"def htheis(r, t, T, S, Q):\n",
" return -Q / (4 * np.pi * T) * exp1(S * r ** 2 / (4 * T * t))\n",
"\n",
"def Qrtheis(r, t, T, S):\n",
" return -Q / (2 * np.pi) * np.exp(-S * r ** 2 / (4 * T * t)) / r"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# solution\n",
"from scipy.special import k0, k1\n",
"def hbar_wbs_res(p, r, T, S, Q, rw, cw, H): # hbar\n",
" sqrtw = np.sqrt(S * p / T)\n",
" Aw = np.pi * rw ** 2\n",
" B = -Q / p * 1 / (p * Aw * k0(rw * sqrtw) + \n",
" (2 * np.pi * rw + Aw * p * cw / H) * \n",
" T * sqrtw * k1(rw * sqrtw))\n",
" hbar = B * k0(r * sqrtw)\n",
" return hbar\n",
"\n",
"def Qrbar_wbs_res(p, r, T, S, Q, rw, cw, H): # Qrbar\n",
" sqrtw = np.sqrt(S * p / T)\n",
" Aw = np.pi * rw ** 2\n",
" B = -Q / p * 1 / (p * Aw * k0(rw * sqrtw) + \n",
" (2 * np.pi * rw + Aw * p * cw / H) * \n",
" T * sqrtw * k1(rw * sqrtw))\n",
" Qrbar = T * sqrtw * B * k1(r * sqrtw)\n",
" return Qrbar\n",
"\n",
"def hwbar_wbs_res(p, r, T, S, Q, rw, cw, H): # hbar inside well\n",
" # r is not used, but Stehfest routine requires it\n",
" hrw = hbar_wbs_res(p, rw, T, S, Q, rw, cw, H)\n",
" Qrw = Qrbar_wbs_res(p, rw, T, S, Q, rw, cw, H)\n",
" return cw * Qrw / H + hrw\n",
"\n",
"# import stehfest_coef and stehfest from Section 5.4\n",
"t = np.logspace(-3, 0, 100)\n",
"ht = htheis(rw, t, T, S, Q) # function from Section 9.1\n",
"hw = stehfest(0, t, hwbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=cw, H=H)\n",
"hwbs = stehfest(rw, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=0, H=H)\n",
"hwbsres = stehfest(rw, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=cw, H=H)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"hideCode": false,
"hidePrompt": false,
"tags": [
"remove-output"
]
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# plot of head at well\n",
"plt.semilogx(t, ht, label='x wbs, x skin')\n",
"plt.semilogx(t, hwbs[0], label=f'v wbs, x sin')\n",
"plt.semilogx(t, hwbsres[0], label=f'v wbs, v skin')\n",
"plt.semilogx(t, hw[0], '--', label=f'in well, v wbs, v skin')\n",
"plt.legend();"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"hideCode": false,
"hidePrompt": false,
"tags": [
"remove-input",
"remove-output"
]
},
"outputs": [
{
"data": {
"application/papermill.record/image/svg+xml": "\n\n\n",
"application/papermill.record/text/plain": ""
},
"metadata": {
"scrapbook": {
"mime_prefix": "application/papermill.record/",
"name": "fig_wellbore_storage"
}
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=diagrams_2d.figdef2)\n",
"plt.subplot(121)\n",
"t = np.logspace(-3, 0, 100)\n",
"ht = htheis(rw, t, T, S, Q)\n",
"hw = stehfest(0, t, hwbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=cw, H=H)\n",
"hwbs = stehfest(rw, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=0, H=H)\n",
"hres = stehfest(rw, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=cw, H=H)\n",
"plt.semilogx(t, ht, label='$r=r_w$, ✘ wbs, ✘ skin')\n",
"plt.semilogx(t, hwbs[0], label=f'$r=r_w$, ✔ wbs, ✘ skin')\n",
"plt.semilogx(t, hres[0], label=f'$r=r_w$, ✔ wbs, ✔ skin')\n",
"plt.semilogx(t, hw[0], '--', label=f'in well, ✔ wbs, ✔ skin')\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.xlabel('time (d)')\n",
"plt.ylabel('head (m)')\n",
"plt.title('(a)', loc='left')\n",
"plt.ylim(-6, 0)\n",
"plt.subplot(122)\n",
"r = 10 # location\n",
"t = np.logspace(-3, 0, 100)\n",
"ht = htheis(r, t, T, S, Q)\n",
"hwbs = stehfest(r, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=0, H=H)\n",
"hres = stehfest(r, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, cw=cw, H=H)\n",
"plt.semilogx(t, ht, label='$r=10$, ✘ wbs, ✘ skin')\n",
"plt.semilogx(t, hwbs[0], label='$r=10$, ✔ wbs, ✘ skin')\n",
"plt.semilogx(t, hres[0], label='$r=10$, ✔ wbs, ✔ skin')\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.xlabel('time (d)')\n",
"plt.ylabel('head (m)')\n",
"plt.title('(b)', loc='left')\n",
"plt.tight_layout()\n",
"glue(\"fig_wellbore_storage\", fig, display=False)"
]
},
{
"cell_type": "markdown",
"metadata": {
"hideCode": false,
"hidePrompt": false
},
"source": [
"```{glue:figure} fig_wellbore_storage\n",
"---\n",
"name: fig_wellbore_storage9.4\n",
"---\n",
"The effect of wellbore storage (wbs) and skin effect (skin) during the first day of pumping; head at the well (left) and head at $r = 10$ m (right).\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"hideCode": false,
"hidePrompt": false
},
"source": [
"```{exercise}\n",
":label: ex_99_head_well\n",
"\n",
"Analyze the pumping test at Oude Korendijk ({numref}`content:chap9.1:theis_well`) taking into account wellbore storage using a well radius of 0.3 m, but no skin effect. Compare your answer to the results of {numref}`content:chap9.1:theis_well`.\n",
"\n",
"\n",
" Solution
\n",
"\n",
"Theis: $T = 560$ m$^2$/d, $S = 4.16 \\cdot 10^{-5}$\n",
"\n",
" \n",
"\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"# # solution\n",
"# from scipy.special import k0, k1\n",
"# def hbar_wbs_res(p, r, T, S, Q, rw, rc, cw): # hbar\n",
"# sqrtw = np.sqrt(S * p / T)\n",
"# A = -Q / p * 1 / (p * np.pi * rc ** 2 * k0(rw * sqrtw) + \n",
"# (2 * np.pi * rw + np.pi * rc ** 2 * p * cw / H) * \n",
"# T * sqrtw * k1(rw * sqrtw))\n",
"# hbar = A * k0(r * sqrtw)\n",
"# return hbar\n",
"\n",
"# def Qrbar_wbs_res(p, r, T, S, Q, rw, rc, cw): # Qrbar\n",
"# sqrtw = np.sqrt(S * p / T)\n",
"# A = -Q / p * 1 / (p * np.pi * rc ** 2 * k0(rw * sqrtw) + \n",
"# (2 * np.pi * rw + np.pi * rc ** 2 * p * cw / H) * \n",
"# T * sqrtw * k1(rw * sqrtw))\n",
"# Qrbar = T * sqrtw * A * k1(r * sqrtw)\n",
"# return Qrbar\n",
"\n",
"# def hwbar_wbs_res(p, r, T, S, Q, rw, rc, cw): # hbar inside well\n",
"# # r is not used, but Stehfest routine requires it\n",
"# hrw = hbar_wbs_res(p, rw, T, S, Q, rw, rc, cw)\n",
"# Qrw = Qrbar_wbs_res(p, rw, T, S, Q, rw, rc, cw)\n",
"# return cw * Qrw / H + hrw\n",
"\n",
"# t = np.logspace(-3, 0, 100)\n",
"# ht = htheis(rw, t, T, S, Q) # function from Section 9.1\n",
"# hw = stehfest(0, t, hwbar_wbs_res, T=T, S=S, Q=Q, rw=rw, rc=rc, cw=cw)\n",
"# hwbs = stehfest(rw, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, rc=rc, cw=0)\n",
"# hwbsres = stehfest(rw, t, hbar_wbs_res, T=T, S=S, Q=Q, rw=rw, rc=rc, cw=cw)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"hideCode": false,
"hidePrompt": false,
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"## Test using ttim\n",
"# k = T / H\n",
"# Ss = S / H\n",
"# import ttim\n",
"# ml = ttim.ModelMaq(kaq=[k], z=[H, 0], Saq=S / H, tmin=1e-3, tmax=10)\n",
"# w = ttim.Well(ml, 0, 0, rw=rw, tsandQ=[(0, Q)], rc=rc, res=cw)\n",
"# ml.solve()\n",
"# httim = ml.head(rw, 0, t)\n",
"# hwttim = w.headinside(t)\n",
"# plt.semilogx(t, httim[0])\n",
"# plt.semilogx(t, hwbsres[0], '--')\n",
"# plt.semilogx(t, hwttim[0])\n",
"# plt.semilogx(t, hw[0], '--');"
]
}
],
"metadata": {
"celltoolbar": "Tags",
"hide_code_all_hidden": false,
"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.13.7"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": false,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}