{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [ "remove-input" ] }, "outputs": [], "source": [ "# %matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\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\n", "import matplotlib as mpl\n", "from matplotlib.lines import Line2D\n", "mpl.rcParams['contour.negative_linestyle'] = 'solid'\n", "from myst_nb import glue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(content:chap10:interface_flow_coast)=\n", "# Interface flow towards the coast" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Several solutions for interface flow in coastal aquifers were discussed in this book so far ({numref}`Chapter %s ` and {numref}`content:chap7:well_near_coast`), all based on the Dupuit approximation. The width of the outflow zone at the coast was zero for all presented Dupuit solutions, except for the solution with a resistance layer along the sea bottom ({numref}`content:chap4:semiconfined_interface_flow`). In {numref}`Chapter %s `, it was claimed that the resulting interface was, nevertheless, a good approximation, even though flow near the coast has a significant vertical component. In this and the following section, Dupuit solutions are compared to solutions for two-dimensional flow in the vertical plane, and it is investigated how wide the outflow zone below the sea is.\n", "
\n", "Consider steady two-dimensional interface flow in the vertical plane in a confined coastal aquifer. The aquifer is so deep that the interface does not reach the bottom of the aquifer; the saltwater is at rest ({numref}`Figure %s `). The density of the freshwater is $\\rho_f$, the density of the saltwater is $\\rho_s$, and the hydraulic conductivity of the aquifer is isotropic and equal to $k$. The sea bottom is horizontal and at the same elevation as the top of the confined aquifer. The origin of the coordinate system is chosen at the coastline at the top of the aquifer. Flow towards the coast is constant and equal to $Q_x=U$. The stream function along the top of the aquifer is chosen to equal $\\psi=0$, so that the stream function along the interface equals $\\psi=U$. \n", "
\n", "This problem is known as a free boundary problem. The location of the interface is a free boundary because it is unknown prior to solving the problem. There are two conditions specified along the interface: the interface is an impermeable boundary and the elevation of the interface $z_i$ is related to the head $h_i$ along the interface according to the Ghyben-Herzberg equation ({numref}`Chapter %s `)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "z_i = -\\alpha h_i\n", "$$(eq_zi10.5)\n", "\n", "where \n", "\n", "$$\n", "\\alpha = \\frac{\\rho_f}{\\rho_s - \\rho_f}\n", "$$(eq_GHconditions10.5)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "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_glover" } }, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = diagrams_2d.glover()\n", "glue(\"fig_glover\", fig, display=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{glue:figure} fig_glover\n", "---\n", "name: fig_glover10.5\n", "---\n", "Interface flow towards the coast.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
The stated problem is known as Glover's problem ({cite:alp}`glover1959pattern`). The complex potential solution used here is obtained from Bruggeman ({cite:author}`bruggeman1999analytical` {cite:year}`bruggeman1999analytical`, Eq. 910.01), as\n", "\n", "$$\n", "\\omega = \\sqrt{-2 k U \\zeta / \\alpha}\n", "$$(eq_omegaglover10.5)\n", "\n", "where $\\zeta=x + iz$ is the complex coordinate. \n", "The solution for the complex potential can be inverted, which gives\n", "\n", "$$\n", "\\zeta = \\frac{-\\omega^2 \\alpha}{2k U}\n", "$$(eq_zetaomegaglover10.5)\n", "\n", "The specific discharge potential along the interface equals $\\phi_i=-k z_i / \\alpha$ and the stream function is equal to the stream function along the bottom of the aquifer ($\\psi=U$). An equation for the elevation of the interface as a function of $x$ may be obtained by substituting $\\omega=-k\\alpha z_i + Ui$ for $\\omega$ and $x+iz_i$ for $\\zeta$ in ({eq}`eq_omegaglover10.5`), which gives, after some algebra\n", "\n", "$$\n", "z_i = -\\sqrt{\\frac{U^2\\alpha^2}{k^2} - \\frac{2U\\alpha x}{k}}\n", "$$(eq_interfaceglover10.5)\n", "\n", "The tip of the interface is the point where the interface intersects the sea bottom ({numref}`Figure %s `). The location of the tip of the interface $x_\\text{tip}$ is obtained by setting $z_i=0$ in Eq. ({eq}`eq_interfaceglover10.5`) and solving for $x$, which gives\n", "\n", "$$\n", "x_\\text{tip} = \\frac{U\\alpha}{2k}\n", "$$(eq_xtip10.5)\n", "\n", "This means that the outflow zone is wider for a larger value of the flow towards the coast $U$. \n", "
In the example below, a flow net is computed by specifically computing the locations of the streamlines and equipotentials, as $\\zeta$ is known as a function of $\\omega$ (Eq. {eq}`eq_zetaomegaglover10.5`) . For the chosen parameters, the interface is already 7 m deep at a distance of 14 m inland from the coastline ({numref}`Figure %s `). The length of the outflow face is just 0.8 m." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# parameters\n", "k = 10 # hydraulic conductivity, m/d\n", "rhof = 1000 # density freshwater, kg/m^3\n", "rhos = 1025 # density of saltwater, kg/m^3\n", "U = 0.4 # flow towards the coast, m^2/d" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "location of the tip: 0.80 m\n" ] } ], "source": [ "# solution\n", "alpha = rhof / (rhos - rhof)\n", "n = 100\n", "zeta_streamlines = []\n", "zeta_equipotentials = []\n", "for psi in np.linspace(0, U, 6):\n", " omega = np.linspace(0, 2, n) + 1j * psi\n", " zeta = -omega**2 * alpha / (2 * k * U)\n", " zeta_streamlines.append(zeta)\n", "for phi in np.arange(0, 2, U / 6):\n", " omega = phi + 1j * np.linspace(0, U, n)\n", " zeta = -omega**2 * alpha / (2 * k * U)\n", " zeta_equipotentials.append(zeta)\n", "zetatip = -(1j * U) ** 2 * alpha / (2 * k * U)\n", "print(f'location of the tip: {zetatip.real: .2f} m')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [ "remove-output" ] }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# basic plot\n", "plt.subplot(111, aspect=1, xlim=(-15, 1))\n", "for zeta in zeta_equipotentials:\n", " plt.plot(zeta.real, zeta.imag, 'C0')\n", "for zeta in zeta_streamlines:\n", " plt.plot(zeta.real, zeta.imag, 'C1')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "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_grovsolfreshwater" } }, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot\n", "fig = plt.figure(figsize=diagrams_2d.figdef2)\n", "plt.subplot(111, aspect=1, xlim=(-15, 1))\n", "for zeta in zeta_equipotentials:\n", " plt.plot(zeta.real, zeta.imag, 'C0')\n", "for zeta in zeta_streamlines[1:-1]:\n", " plt.plot(zeta.real, zeta.imag, 'C1')\n", "zeta = zeta_streamlines[0]\n", "plt.plot(zeta.real, zeta.imag, 'k') \n", "zeta = zeta_streamlines[-1]\n", "plt.plot(zeta.real, zeta.imag, color=diagrams_2d.C1dark)\n", "diagrams_2d.arrow((0, 0.2), (zetatip.real, 0.2), arrow='<->')\n", "diagrams_2d.arrow((0, -3), (0, -5), 'coast \\n line', arrow='->', ha='center')\n", "plt.text(0.4, 0.4, 'outflow zone', ha='center')\n", "plt.ylim(-7.5, 1)\n", "plt.xlim(-15, 2.5)\n", "plt.xlabel('x (m)')\n", "plt.ylabel('z (m)')\n", "plt.tight_layout()\n", "glue('fig_grovsolfreshwater',fig ,display=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{glue:figure} fig_grovsolfreshwater\n", "---\n", "name: fig_grovsolfreshwater10.5\n", "---\n", "Flow net for two-dimensional flow of fresh water towards the coast (Glover's solution).\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "     \n", "The interface position obtained with Glover's solution is compared to a Dupuit interface solution. The Dupuit solution for the elevation of the interface is obtained from {numref}`content:chap4:confined_flow` as\n", "\n", "$$\n", "z_i = -\\sqrt{\\frac{-2U \\alpha x}{k}}\n", "$$(eq_zisqrtfrac10.5)\n", "\n", "Note that $z_i^2$ of the two-dimensional solution ({eq}`eq_interfaceglover10.5`) differs by a constant factor from the Dupuit solution. The difference between the Glover solution and the Dupuit solution is minor for the parameters of the example ({numref}`Figure %s `). \n", "
      The difference between the two solutions becomes a bit larger when the vertical hydraulic conductivity $k_z$ of the aquifer is significantly smaller than the horizontal hydraulic conductivity $k_x$. The Glover solution is modified for an anistropic hydraulic conductivity using the approach of {numref}`content:chap10_vertical_anisotropy`. The solution for the complex potential is \n", "\n", "$$\n", "\\tilde{\\omega} = \\sqrt{-2 \\tilde{k} U \\tilde{\\zeta} / \\alpha}\n", "$$(eq_omegatildesqrt10.5)\n", "\n", "The corresponding elevation of the interface is\n", "\n", "$$\n", "z_i = -\\sqrt{\\frac{U^2\\alpha^2}{k_x k_z} - \\frac{2U\\alpha x}{k_x}}\n", "$$(eq_zialpha10.5)\n", "\n", "and the location of the tip of the interface is\n", "\n", "$$\n", "x_\\text{tip} = \\frac{U\\alpha}{2k_z}\n", "$$(eq_xtipalpha10.5)\n", "\n", "In the example below, the interface is plotted for the isotropic case and for the case that $k_z=k_x/10$, which results in an outflow zone that is 10 times wider than for the isotropic case ({numref}`Figure %s `). The interface is located somewhat deeper than for the isotropic case. The Dupuit solution does not change, of course, as the resistance to flow in the vertical direction is neglected. \n", "
      The Dupuit solution is a very good approximation of the isotropic Glover solution and a reasonable approximation of the anisotropic Glover solution. The Dupuit solution is always conservative when assessing the extent of seawater intrusion: the thickness of the freshwater zone is always smaller and the seawater intrusion more severe than in the Glover solution. The largest difference in the thickness of the freshwater zone is at the coastline." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# additional parameters anisotropic case\n", "kx = k\n", "kz = k / 10" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "location of the tip for anisotropic case: 8.00 m\n" ] } ], "source": [ "# isotropic solution\n", "xtip = U * alpha / (2 * k)\n", "xglover = np.linspace(-150, xtip, 100)\n", "zglover = -np.sqrt(U * alpha / k * (U * alpha / k - 2 * xglover))\n", "# aniosotropic solution\n", "xtipaniso = U * alpha / (2 * kz)\n", "xaniso = np.linspace(-150, xtipaniso, 200)\n", "zaniso = -np.sqrt(U * alpha / kx * (U * alpha / kz - 2 * xaniso))\n", "print(f'location of the tip for anisotropic case: {xtipaniso: 0.2f} m')\n", "#\n", "xdupuit = np.linspace(-150, 0, 100)\n", "zdupuit = -np.sqrt(-2 * U * alpha * xdupuit / k)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [ "remove-output" ] }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# basic plot (Fig. 10.15)\n", "plt.plot(xglover, zglover, 'C1', label='Glover $k_z=k_x$')\n", "plt.plot(xaniso, zaniso, 'C2', label='Glover $k_z=k_x/10$')\n", "plt.plot(xdupuit, zdupuit, 'C1--', label='Dupuit $H=20$')\n", "plt.legend();" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [ "remove-output", "remove-input" ] }, "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_grovdupuit" } }, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot\n", "fig = plt.figure(figsize=diagrams_2d.figdef1new)\n", "plt.subplot(111, aspect=4, xlim=(-150, 10), ylim=(-24, 1))\n", "plt.plot(xglover, zglover, diagrams_2d.C1dark, label='Glover $k_z=k_x$')\n", "plt.plot(xaniso, zaniso, 'C1', label='Glover $k_z=k_x/10$')\n", "plt.plot(xdupuit, zdupuit, '--', color=diagrams_2d.C1dark, label='Dupuit')\n", "plt.xlabel('x (m)')\n", "plt.ylabel('z (m) - VE:4x')\n", "plt.legend()\n", "plt.grid()\n", "plt.tight_layout();\n", "glue('fig_grovdupuit',fig ,display=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{glue:figure} fig_grovdupuit\n", "---\n", "name: fig_grovdupuit10.5\n", "---\n", "Comparison of interface elevation between the Glover solution (both isotropic and anisotropic) and the Dupuit solution.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{exercise}\n", ":label: ex_10.7\n", "Compute the ratio of the depth of the interface and the head at the interface at $x=-50$ and at $x=-100$ for the anisotropic case (Hint: the answer is of course $\\alpha$ if you do this correctly). Next, compute the ratio of the depth of the interface and the head at the top of the aquifer at $x=-50$ and at $x=-100$ for the anisotropic case. \n", "\n", "

\n", "

Solution

\n", "\n", "Factor based on head at the top of aquifer at $x = −100$ m: $41.6$. \\\n", "Factor based on head at the top of aquifer at $x =−50$ m: $43.1$\n", "\n", "

\n", "\n", "```\n" ] } ], "metadata": { "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 } }, "nbformat": 4, "nbformat_minor": 4 }