{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"# %matplotlib inline\n",
"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 matplotlib as mpl\n",
"mpl.rcParams['contour.negative_linestyle'] = 'solid'"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"(content:chap10_vertical_anisotropy)=\n",
"# Vertical anisotropy\n",
"\n",
"The hydraulic conductivity is often anisotropic, where the vertical component of the hydraulic conductivity is smaller than the horizontal component. One of the main reasons for vertical anisotropy is that aquifers often consist of alternating layers of somewhat coarser and somewhat finer material. The average value of the hydraulic conductivity in the horizontal direction is dictated by the hydraulic conductivity of the coarser layers, while the average value in the vertical direction is dictated by the hydraulic conductivity of the finer layers. In practice, it is popular to use an anisotropy factor of 10 (the vertical hydraulic conductivity is 10 times smaller than the horizontal hydraulic conductivity), but a good reason for this number does not exist. For the case that the hydraulic conductivity is anisotropic with $k_x$ in the $x$-direction and $k_z$ in the $z$-direction, Darcy's law becomes \n",
"\n",
"\n",
"$$\n",
"q_x = -k_x\\frac{\\partial h}{\\partial x} {\\hskip 2em}\n",
"q_z = -k_z\\frac{\\partial h}{\\partial z}\n",
"$$(eq_darcyslaw10.1)\n",
"\n",
"Substitution of these equations for $q_x$ and $q_z$ in the continuity equation ({eq}`eq_2dsteadyverticalcont`) gives\n",
"\n",
"$$\n",
"\\frac{\\partial}{\\partial x}\\left(k_x \\frac{\\partial h}{\\partial x}\\right) + \n",
"\\frac{\\partial}{\\partial z}\\left(k_z \\frac{\\partial h}{\\partial z}\\right)= 0\n",
"$$(eq_qxqz10.1)\n",
"\n",
"This equation may be transformed into Laplace's equation by a change of variables\n",
"\n",
"$$\n",
"\\tilde{x} = x\\sqrt{k_z/k_x} {\\hskip 2em} \\tilde{z} = z\n",
"$$(eq_anisosqueeze10.1)\n",
"\n",
"which gives\n",
"\n",
"$$\n",
"\\frac{\\partial^2 h}{\\partial \\tilde{x}^2} + \\frac{\\partial^2 h}{\\partial \\tilde{z}^2} = 0\n",
"$$(eq_partial10.1)\n",
"\n",
"Hence, in the $\\tilde{x}, \\tilde{z}$ domain, the head is governed by Laplace's equation. The specific discharge in the transformed domain can be evaluated with Darcy's law using the isotropic hydraulic conductivity $\\tilde{k}$ defined as\n",
"\n",
"$$\n",
"\\tilde{k} = \\sqrt{k_xk_z}\n",
"$$(eq_isok10.1)\n",
"\n",
"The use of $\\tilde{k}$ ensures that the net discharge in the isotropic transformed domain is the same as the net discharge in the original anisotropic domain. For example, the net flux $Q_\\text{net}$ across a small section $\\Delta z$ in the anisotropic domain equals\n",
"\n",
"$$\n",
"Q_\\text{net} = -k_x \\frac{\\partial h}{\\partial x} \\Delta z\n",
"$$(eq_Qnet10.1)\n",
"\n",
"and is equal to the net flux across the same section in the transformed isotropic domain\n",
"\n",
"$$\n",
"\\tilde{Q}_\\text{net} = -\\tilde{k}\\frac{\\partial h}{\\partial \\tilde{x}} \\Delta \\tilde{z}=\n",
"-\\sqrt{k_xk_z}\\frac{\\partial h}{\\partial x} \\frac{1}{\\sqrt{k_z/k_x}} \\Delta z = -k_x \\frac{\\partial h}{\\partial x} \\Delta z = Q_\\text{net}\n",
"$$(eq_Qnettilde10.1)\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"The same holds in the $z$-direction.\n",
" \n",
"In summary, an anisotropic problem may be solved with standard methods for Laplace's equation, provided that the $x$ direction of the flow field is squeezed according to Eq. ({eq}`eq_anisosqueeze10.1`), and the isotropic hydraulic conductivity (Eq. {eq}`eq_isok10.1`) is used in the transformed domain. A complex specific discharge potential in the transformed potential may be defined as $\\tilde{\\omega}=\\tilde{\\phi} + i\\tilde{\\psi}$, where $\\phi=\\tilde{k}h$ and $\\tilde{\\psi} = \\psi$. \n",
" \n",
"As a simple example, consider a rectangular domain with length $L$ and height $B$ ({numref}`Figure %s `). The head at the origin in the lower left-hand corner equals zero. Flow is uniform, and the head gradients in the $x$- and $z$-directions are equal: $\\partial h/\\partial x= \\partial h / \\partial z = a$. As such, the head varies linearly along the boundaries ({numref}`Figure %s `). The hydraulic conductivities in the $x$ and $z$ directions are $k_x$ and $k_z$, respectively. \n",
" \n",
"The flow is uniform in both the anisotropic physical domain and the isotropic transformed domain. In the physical domain, $q_x=-k_xa$ and $q_z=-k_za$. In the transformed domain, $\\tilde{q}_x$ and $\\tilde{q}_z$ are\n",
"\n",
"$$\n",
"\\tilde{q}_x = -\\tilde{k}\\frac{\\partial h}{\\partial \\tilde{x}} = -\\tilde{k}\\frac{\\partial h}{\\partial x}\\frac{\\partial x}{\\partial \\tilde{x}} = -\\frac{\\tilde{k}a}{\\sqrt{k_z/k_x}}\n",
"$$(eq_qxtilde10.1)\n",
"\n",
"$$\n",
"\\tilde{q}_z=-\\tilde{k}a\n",
"$$(eq_qztilde10.1)\n",
"\n",
"The solution for the complex specific discharge potential in the transformed domain is (see {numref}`content:chap7:well_near_river`)\n",
"\n",
"$$\n",
"\\tilde{\\omega} = -(\\tilde{q}_x- \\tilde{q}_zi) \\tilde{\\zeta}\n",
"$$(eq_omegatilde10.1)\n",
"\n",
"where $\\tilde{\\zeta}=\\tilde{x} + i\\tilde{z}$ is the complex coordinate in the transformed domain. \n",
" \n",
"An example flow field is shown below in {numref}`Figure %s `. The flow field is shown in both the transformed domain and in the problem domain. Note that the streamlines and head contours intersect at right angles in the transformed domain, where the hydraulic conductivity $\\tilde{k}$ is isotropic, but don't intersect at right angles in the problem domain where the hydraulic conductivity is anisotropic. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"# parameters\n",
"L = 1000 # length of domain, m\n",
"B = 300 # width of domain, m\n",
"kx = 5 # hydraulic conductivity in x-direction, m/d\n",
"kz = 1 # hydraulic conductivity in z-direction, m/d\n",
"a = -0.002 # head gradient in x and z directions"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"# solution\n",
"kt = np.sqrt(kx * kz)\n",
"qxt = -kt * a / np.sqrt(kz / kx)\n",
"qzt = -kt * a\n",
"xg, zg = np.meshgrid(np.linspace(0, L, 10), np.linspace(0, B, 10))\n",
"xt = xg * np.sqrt(kz / kx)\n",
"zt = zg\n",
"zetat = xt + 1j * zt\n",
"omt = -(qxt - 1j * qzt) * zetat\n",
"h = omt.real / kt\n",
"psi = omt.imag"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"remove-output"
]
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# basic flow net\n",
"plt.subplot(121, aspect=1) # problem domain\n",
"plt.contour(xg, zg, h, 10)\n",
"plt.contour(xg, zg, psi, 10);\n",
"plt.subplot(122, aspect=1) # transformed domain\n",
"plt.contour(xt, zt, h, 20)\n",
"plt.contour(xt, zt, psi, 20);"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"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_uniflowanisotropic"
}
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from matplotlib import gridspec\n",
"from matplotlib.lines import Line2D\n",
"# Create figure and grid spec\n",
"fig = plt.figure(figsize=(8, 4))\n",
"spec = gridspec.GridSpec(ncols=2, nrows=2, width_ratios=[2.25, 1])\n",
"\n",
"# --- Plot a ---\n",
"ax0 = fig.add_subplot(spec[0], aspect=1)\n",
"ax0.set_title('(a)', loc='left')\n",
"ax0.plot([0, L, L, 0, 0], [0, 0, B, B, 0], 'k')\n",
"ax0.text(L/2, 10, '$h=ax$', ha='center', va='bottom')\n",
"ax0.text(L/2, B-10, '$h=ax + aB$', ha='center', va='top')\n",
"ax0.text(10, B/2, '$h=az$', rotation=90, va='center')\n",
"ax0.text(L-10, B/2, '$h=az + aL$', rotation=90, va='center', ha='right')\n",
"ax0.set_xlim(0, L)\n",
"ax0.set_ylim(0, B)\n",
"ax0.set_xticks([0, L])\n",
"ax0.set_xticklabels(['0', '$L$'])\n",
"ax0.set_yticks([0, B])\n",
"ax0.set_yticklabels(['0', '$B$'])\n",
"ax0.set_xlabel('$x$')\n",
"ax0.set_ylabel('$z$')\n",
"\n",
"# Define contour levels\n",
"h_levels = np.arange(-2.4, 0, 0.4)\n",
"psi_levels = np.arange(-6, 2, 1)\n",
"\n",
"# --- Plot b ---\n",
"ax1 = fig.add_subplot(spec[2], aspect=1)\n",
"ax1.set_title('(b)', loc='left')\n",
"\n",
"cs_h = ax1.contour(xg, zg, h, h_levels, colors='C0')\n",
"ax1.clabel(cs_h, fmt='%1.1f')\n",
"\n",
"cs_psi = ax1.contour(xg, zg, psi, psi_levels, colors='C1')\n",
"ax1.clabel(cs_psi, fmt='%1.1f')\n",
"\n",
"ax1.set_xlabel('$x$')\n",
"ax1.set_ylabel('$z$')\n",
"\n",
"# Add legend via proxy artists\n",
"legend_elements = [\n",
" Line2D([0], [0], color='C0', label='$h$'),\n",
" Line2D([0], [0], color='C1', label=r'$\\psi$')\n",
"]\n",
"ax1.legend(handles=legend_elements, loc='upper right')\n",
"\n",
"# --- Plot c ---\n",
"ax2 = fig.add_subplot(spec[3], aspect=1)\n",
"ax2.set_title('(c)', loc='left')\n",
"\n",
"cs_h_t = ax2.contour(xt, zt, h, h_levels, colors='C0')\n",
"ax2.clabel(cs_h_t, fmt='%1.1f')\n",
"\n",
"cs_psi_t = ax2.contour(xt, zt, psi, psi_levels, colors='C1')\n",
"ax2.clabel(cs_psi_t, fmt='%1.1f', fontsize='smaller')\n",
"\n",
"ax2.set_xlabel(r'$\\tilde{x}$')\n",
"ax2.set_ylabel(r'$\\tilde{z}$')\n",
"\n",
"# Same legend pattern for this subplot\n",
"ax2.legend(handles=legend_elements, loc='upper right')\n",
"\n",
"plt.tight_layout()\n",
"glue('fig_uniflowanisotropic', fig, display=False)\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"```{glue:figure} fig_uniflowanisotropic\n",
"---\n",
"name: fig_aniuniflow10.1\n",
"---\n",
"Example for uniform flow in an anisotropic aquifer; boundary conditions (upper left), flow net in squeezed transformed domain (lower right), and flow net in problem domain (lower left).\n",
"```"
]
}
],
"metadata": {
"author": "k",
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "jupbook",
"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.2"
},
"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
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}