{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Building Pit\n", "\n", "This example demonstrates the modeling of a building pit with dewatering wells using `timflow`. The model calculates groundwater flow towards the pit and evaluates the effectiveness of dewatering strategies by computing total discharge and drawdown distances.\n", "\n", "
\n", "\n", "Note\n", "\n", "It is highly recommended to use the BuildingPit element if you want to\n", "implement different layer properties inside the building pit as compared to the rest of\n", "the aquifer. Adding barriers, e.g. `LeakyWalls`, around an inhomogeneity will\n", "not work.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "import timflow.steady as tfs" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Parameters\n", "\n", "Define some aquifer parameters " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kh = 2.0 # m/day\n", "f_ani = 1 / 10 # anisotropy factor\n", "kv = f_ani * kh\n", "ctop = 800.0 # resistance top leaky layer in days" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define aquifer top and bottom, the depth of the sheetpile wall and the position of the dewatering wells." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ztop = 0.0 # surface elevation\n", "z_well = -13.0 # end depth of the wellscreen\n", "z_dw = -15.0 # bottom elevation of sheetpile wall\n", "z_extra = z_dw - 15.0 # extra layer\n", "zbot = -60.0 # bottom elevation of the model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Size of the building pit, and the required drawdown at the center of the building pit." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "length = 40.0 # length building pit in m\n", "width = 30.0 # width building pit in m\n", "\n", "h_bem = -6.21 # m\n", "offset = 5.0 # distance groundwater extraction element from sheetpiles in m" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xy = [\n", " (-length / 2, -width / 2),\n", " (length / 2, -width / 2),\n", " (length / 2, width / 2),\n", " (-length / 2, width / 2),\n", " (-length / 2, -width / 2),\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the building pit" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "(p2,) = plt.plot(*np.array(xy).T, \"o-\")\n", "plt.axis(\"equal\")\n", "plt.show()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Model\n", "Set up a model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "z = np.array([ztop + 1, ztop, z_dw, z_dw, z_extra, z_extra, zbot])\n", "dz = z[1::2] - z[2::2]\n", "kh_arr = kh * np.ones(dz.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Resistances of the top confining layer and aquitards" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c = np.r_[np.array([ctop]), dz[:-1] / (2 * kv) + dz[1:] / (2 * kv)]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ml = tfs.ModelMaq(kaq=kh_arr, z=z, c=c, topboundary=\"semi\", hstar=0.0)\n", "\n", "layers = np.arange(np.sum(z_dw <= ml.aq.zaqbot))\n", "last_lay_dw = layers[-1]\n", "\n", "inhom = tfs.BuildingPitMaq(\n", " ml,\n", " xy,\n", " kaq=kh_arr,\n", " z=z[1:],\n", " topboundary=\"conf\",\n", " c=c[1:],\n", " order=4,\n", " ndeg=3,\n", " layers=layers,\n", ")\n", "\n", "tfs.River(\n", " ml,\n", " x1=-length / 2 + offset,\n", " y1=width / 2 - offset,\n", " x2=length / 2 - offset,\n", " y2=width / 2 - offset,\n", " hls=h_bem,\n", " layers=np.arange(last_lay_dw + 1),\n", ")\n", "tfs.River(\n", " ml,\n", " x1=-length / 2 + offset,\n", " y1=0,\n", " x2=length / 2 - offset,\n", " y2=0,\n", " hls=h_bem,\n", " layers=np.arange(last_lay_dw + 1),\n", ")\n", "tfs.River(\n", " ml,\n", " x1=-length / 2 + offset,\n", " y1=-width / 2 + offset,\n", " x2=length / 2 - offset,\n", " y2=-width / 2 + offset,\n", " hls=h_bem,\n", " layers=np.arange(last_lay_dw + 1),\n", ")\n", "\n", "# ml.solve_mp(nproc=2)\n", "ml.solve()\n", "\n", "Qtot = 0.0\n", "\n", "for e in ml.elementlist:\n", " if e.name == \"River\":\n", " Qtot += e.discharge()\n", "\n", "print(\"\\nDischarge =\", np.round(Qtot.sum(), 2), \"m^3/day\")\n", "\n", "y = np.linspace(-width / 2 - 25, width / 2 + 1100, 201)\n", "hl = ml.headalongline(np.zeros(201), y, layers=[0])\n", "y_5cm = np.interp(\n", " -0.05, ml.headalongline(np.zeros(201), y, layers=0).squeeze(), y, right=np.nan\n", ")\n", "print(\"Distance to 5 cm drawdown contour =\", np.round(y_5cm, 2), \"m\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Plot an overview of the model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ml.plots.topview()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Visualizations\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(0.0, length / 2 + 1100, 201)\n", "hl = ml.headalongline(x, np.zeros(201), layers=[last_lay_dw, last_lay_dw + 1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n", "ax.plot(x, hl[0].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n", "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n", "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n", "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n", "ax.set_xlabel(\"x (m)\")\n", "ax.set_ylabel(\"head (m)\")\n", "ax.legend(loc=\"best\")\n", "ax.grid(True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(-length / 2 - 25, 0.0, 201)\n", "hl = ml.headalongline(x, np.zeros(201), layers=[last_lay_dw, last_lay_dw + 1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n", "\n", "ax.plot(x, hl[0].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n", "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n", "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n", "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n", "ax.set_xlabel(\"x (m)\")\n", "ax.set_ylabel(\"head (m)\")\n", "ax.legend(loc=\"best\")\n", "ax.grid(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot cross-section around the sheetpile wall" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xoffset = 50\n", "zoffset = 15\n", "x1, x2, y1, y2 = [-length / 2 - xoffset, 0.0, 0, 0]\n", "nudge = 1e-6\n", "n = 301" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot head contour cross-sections\n", "h = ml.headalongline(\n", " np.linspace(x1 + nudge, x2 - nudge, n),\n", " np.linspace(y1 + nudge, y2 - nudge, n),\n", ")\n", "L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)\n", "xg = np.linspace(0, L, n) + x1\n", "\n", "zg = 0.5 * (ml.aq.zaqbot + ml.aq.zaqtop)\n", "zg = np.hstack((ml.aq.zaqtop[0], zg, ml.aq.zaqbot[-1]))\n", "h = np.vstack((h[0], h, h[-1]))\n", "\n", "ymid = np.mean([y1, y2])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "levels = np.linspace(h_bem - 0.1, -0.0, 51)\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n", "ax.set_aspect(\"equal\")\n", "ml.plots.xsection(\n", " xy=[(x1, ymid), (x2, ymid)],\n", " ax=ax,\n", " horizontal_axis=\"x\",\n", " labels=True,\n", " zorder=10,\n", ")\n", "cf = ax.contourf(xg, zg, h, levels)\n", "cs = ax.contour(xg, zg, h, levels, colors=\"k\", linewidths=0.5)\n", "ax.set_ylim(z_dw - zoffset, z_dw + zoffset)\n", "ax.set_ylabel(\"depth (m NAP)\")\n", "ax.set_xlabel(\"m\")\n", "\n", "cb = plt.colorbar(cf, ax=ax, shrink=0.6)\n", "cb.set_label(\"head (m)\")\n", "cb.set_ticks(np.arange(-6, 0.1, 1))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Model 2: Add more layers\n", "Add more layers to the model to get a more accurate solution of the flow towards the building pit." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 11 # number of layers around bottom of sheetpile wall" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculate thickness of each sublayer above and below the sheetpile wall\n", "dz_i_top = (z_well - z_dw) / np.sum(np.arange(n + 1))\n", "dz_i_bot = (z_dw - z_extra) / np.sum(np.arange(2 * n + 1))\n", "\n", "# Generate cumulative depths for sublayers above and below the wall\n", "z_layers_top = np.cumsum(np.arange(n) * dz_i_top)\n", "z_layers_bot = np.cumsum(np.arange(2 * n) * dz_i_bot)\n", "\n", "# Combine sublayer depths into a single array for the model\n", "zgr = np.r_[z_dw + z_layers_top[::-1], (z_dw - z_layers_bot)[1:]]\n", "\n", "# Build full array of layer boundaries for the model\n", "z4 = np.r_[\n", " np.array([ztop + 1, ztop, z_well, z_well]),\n", " np.repeat(zgr, 2, 0),\n", " z_extra,\n", " z_extra,\n", " zbot,\n", "]\n", "\n", "# Calculate thicknesses and hydraulic conductivities for each layer\n", "dz4 = z4[1:-1:2] - z4[2::2]\n", "kh_arr = kh * np.ones(dz4.shape)\n", "\n", "# Calculate resistance for each layer\n", "c = np.r_[np.array([ctop]), dz4[:-1] / (2 * kv) + dz4[1:] / (2 * kv)]\n", "\n", "# Set hydraulic conductivity of the top layer to a very low value\n", "kh_arr2 = kh_arr.copy()\n", "kh_arr2[0] = 1e-5" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ml = tfs.ModelMaq(kaq=kh_arr, z=z4, c=c, topboundary=\"semi\", hstar=0.0)\n", "\n", "layers = np.arange(np.sum(z_dw <= ml.aq.zaqbot))\n", "last_lay_dw = layers[-1]\n", "inhom = tfs.BuildingPitMaq(\n", " ml,\n", " xy,\n", " kaq=kh_arr2,\n", " z=z4[1:],\n", " topboundary=\"conf\",\n", " c=c[1:],\n", " order=4,\n", " ndeg=3,\n", " layers=layers,\n", ")\n", "\n", "wlayers = np.arange(np.sum(-14 <= ml.aq.zaqbot))\n", "wlayers = wlayers[1:]\n", "\n", "tfs.River(\n", " ml,\n", " x1=-length / 2 + offset,\n", " y1=width / 2 - offset,\n", " x2=length / 2 - offset,\n", " y2=width / 2 - offset,\n", " hls=h_bem,\n", " layers=wlayers,\n", ")\n", "tfs.River(\n", " ml,\n", " x1=-length / 2 + offset,\n", " y1=0,\n", " x2=length / 2 - offset,\n", " y2=0,\n", " hls=h_bem,\n", " layers=wlayers,\n", " order=5,\n", ")\n", "tfs.River(\n", " ml,\n", " x1=-length / 2 + offset,\n", " y1=-width / 2 + offset,\n", " x2=length / 2 - offset,\n", " y2=-width / 2 + offset,\n", " hls=h_bem,\n", " layers=wlayers,\n", ")\n", "\n", "# ml.solve_mp(nproc=2)\n", "ml.solve()\n", "\n", "Qtot_ml = 0.0\n", "\n", "for e in ml.elementlist:\n", " if e.name == \"River\":\n", " Qtot_ml += e.discharge()\n", "\n", "print(\"\\nDischarge =\", np.round(Qtot_ml.sum(), 2), \"m^3/day\")\n", "\n", "y = np.linspace(-width / 2 - 25, width / 2 + 1100, 201)\n", "hl = ml.headalongline(np.zeros(201), y, layers=[0])\n", "y_5cm_ml = np.interp(\n", " -0.05, ml.headalongline(np.zeros(201), y, layers=0).squeeze(), y, right=np.nan\n", ")\n", "print(\"Distance to 5 cm drawdown contour =\", np.round(y_5cm_ml, 2), \"m\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "last_lay_dw = layers[-1]\n", "x = np.linspace(0.0, length / 2 + 1100, 201)\n", "hl = ml.headalongline(x, np.zeros(201), layers=[0, last_lay_dw, last_lay_dw + 1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n", "\n", "ax.plot(x, hl[0].squeeze(), label=\"head layer 0\")\n", "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n", "ax.plot(x, hl[2].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n", "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n", "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n", "ax.set_xlabel(\"x (m)\")\n", "ax.set_ylabel(\"head (m)\")\n", "ax.legend(loc=\"best\")\n", "ax.grid(True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(-length / 2 - 25, 0.0, 201)\n", "hl = ml.headalongline(x, np.zeros(201), layers=[0, last_lay_dw, last_lay_dw + 1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n", "\n", "ax.plot(x, hl[0].squeeze(), label=\"head layer 0\")\n", "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n", "ax.plot(x, hl[2].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n", "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n", "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n", "ax.set_xlabel(\"x (m)\")\n", "ax.set_ylabel(\"head (m)\")\n", "ax.legend(loc=\"best\")\n", "ax.grid(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot head contours in a vertical cross-section around the sheetpile wall." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xoffset = 50\n", "zoffset = 15\n", "x1, x2, y1, y2 = [-length / 2 - xoffset, 0.0, 0, 0]\n", "nudge = 1e-6\n", "n = 301" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot head contour cross-sections\n", "h = ml.headalongline(\n", " np.linspace(x1 + nudge, x2 - nudge, n), np.linspace(y1 + nudge, y2 - nudge, n)\n", ")\n", "L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)\n", "xg = np.linspace(0, L, n) + x1\n", "\n", "zg = 0.5 * (ml.aq.zaqbot + ml.aq.zaqtop)\n", "zg = np.hstack((ml.aq.zaqtop[0], zg, ml.aq.zaqbot[-1]))\n", "h = np.vstack((h[0], h, h[-1]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "levels = np.linspace(h_bem - 0.1, -0.0, 51)\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n", "ax.set_aspect(\"equal\")\n", "ml.plots.xsection(xy=[(x1, ymid), (x2, ymid)], ax=ax, labels=False, horizontal_axis=\"x\")\n", "cf = ax.contourf(xg, zg, h, levels)\n", "cs = ax.contour(xg, zg, h, levels, colors=\"k\", linewidths=0.5)\n", "ax.set_ylim(z_dw - zoffset, z_dw + zoffset)\n", "ax.set_ylabel(\"depth (m NAP)\")\n", "ax.set_xlabel(\"m\")\n", "\n", "cb = plt.colorbar(cf, ax=ax, shrink=0.6)\n", "cb.ax.set_ylabel(\"head (m)\")\n", "cb.set_ticks(np.arange(-6, 0.1, 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note the difference between the computed discharge between the model with only a few layers and the model with very fine layers." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Number of layers | N=3 | N=35 | unit\")\n", "print(\"-\" * 56)\n", "print(\n", " f\"Discharge from building pit | {Qtot.sum():>5.1f} \"\n", " f\"| {Qtot_ml.sum():>5.1f} | m^3/d\"\n", ")\n", "print(f\"Distance to 5 cm drawdown contour | {y_5cm:>5.1f} | {y_5cm_ml:>5.1f} | m\")" ] } ], "metadata": { "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 4 }