From cd40298aebe86010338a82c770640ba139b76bfb Mon Sep 17 00:00:00 2001 From: raar1 Date: Wed, 4 Feb 2026 14:30:16 +0100 Subject: [PATCH 01/10] Resolve rebase conflict --- {examples => docs/examples}/additional_term.ipynb | 0 .../examples}/example_equation_formulation.ipynb | 0 {examples => docs/examples}/example_exact_equation.ipynb | 0 .../examples}/example_mibitrans_anatrans_comparison.ipynb | 0 {examples => docs/examples}/example_walkthrough.ipynb | 0 mkdocs.yml | 7 +++++-- pyproject.toml | 1 + 7 files changed, 6 insertions(+), 2 deletions(-) rename {examples => docs/examples}/additional_term.ipynb (100%) rename {examples => docs/examples}/example_equation_formulation.ipynb (100%) rename {examples => docs/examples}/example_exact_equation.ipynb (100%) rename {examples => docs/examples}/example_mibitrans_anatrans_comparison.ipynb (100%) rename {examples => docs/examples}/example_walkthrough.ipynb (100%) diff --git a/examples/additional_term.ipynb b/docs/examples/additional_term.ipynb similarity index 100% rename from examples/additional_term.ipynb rename to docs/examples/additional_term.ipynb diff --git a/examples/example_equation_formulation.ipynb b/docs/examples/example_equation_formulation.ipynb similarity index 100% rename from examples/example_equation_formulation.ipynb rename to docs/examples/example_equation_formulation.ipynb diff --git a/examples/example_exact_equation.ipynb b/docs/examples/example_exact_equation.ipynb similarity index 100% rename from examples/example_exact_equation.ipynb rename to docs/examples/example_exact_equation.ipynb diff --git a/examples/example_mibitrans_anatrans_comparison.ipynb b/docs/examples/example_mibitrans_anatrans_comparison.ipynb similarity index 100% rename from examples/example_mibitrans_anatrans_comparison.ipynb rename to docs/examples/example_mibitrans_anatrans_comparison.ipynb diff --git a/examples/example_walkthrough.ipynb b/docs/examples/example_walkthrough.ipynb similarity index 100% rename from examples/example_walkthrough.ipynb rename to docs/examples/example_walkthrough.ipynb diff --git a/mkdocs.yml b/mkdocs.yml index 3c4f020..6df5e1d 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -11,7 +11,10 @@ nav: - Home: index.md - Introduction: introduction/introduction.md - Background: - - Bioscreen: background/bioscreen.md + - Bioscreen: background/bioscreen.md + - Examples: + - Walkthrough: examples/example_walkthrough.ipynb + - Comparison Mibitrans-Anatrans: examples/example_mibitrans_anatrans_comparison.ipynb - Development: development/development.md - API Reference: - Data: reference/reference_data.md @@ -67,7 +70,7 @@ plugins: show_source: true options: show_submodules: true - +- mkdocs-jupyter # Styled blocks: https://squidfunk.github.io/mkdocs-material/reference/admonitions/#supported-types markdown_extensions: diff --git a/pyproject.toml b/pyproject.toml index 46e0bb8..7ba4311 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -85,6 +85,7 @@ doc = [ "mdx-include >=1.4.1", "mkdocs-markdownextradata-plugin >=0.2.5", "mike >=2.1.3" + "mkdocs-jupyter >=0.25.1" ] [project.urls] From 45d4d1a94a1c75b7a46ce8c117c00a9c2202534d Mon Sep 17 00:00:00 2001 From: raar1 Date: Mon, 19 Jan 2026 13:23:50 +0100 Subject: [PATCH 02/10] Update nb check action --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 11e8de6..e4dd5a7 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -63,4 +63,4 @@ jobs: - name: Check style against standards using ruff run: ruff check . - name: Check cleanliness of notebooks - run: nb-clean check examples/ + run: nb-clean check docs/examples/ From 0cf8949deedb344bc311a4261fd2fe009d18df76 Mon Sep 17 00:00:00 2001 From: raar1 Date: Mon, 19 Jan 2026 15:37:17 +0100 Subject: [PATCH 03/10] Execute notebooks and provide source link --- mkdocs.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/mkdocs.yml b/mkdocs.yml index 6df5e1d..3e91060 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -70,7 +70,9 @@ plugins: show_source: true options: show_submodules: true -- mkdocs-jupyter +- mkdocs-jupyter: + execute: true + include_source: True # Styled blocks: https://squidfunk.github.io/mkdocs-material/reference/admonitions/#supported-types markdown_extensions: From eba92573f94fd9f446fa67aee66382b3cc6579a8 Mon Sep 17 00:00:00 2001 From: raar1 Date: Wed, 4 Feb 2026 14:15:17 +0100 Subject: [PATCH 04/10] Reformat js file --- docs/assets/mathjax.js | 39 ++++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/docs/assets/mathjax.js b/docs/assets/mathjax.js index 35b2081..d5a5cd3 100644 --- a/docs/assets/mathjax.js +++ b/docs/assets/mathjax.js @@ -1,21 +1,22 @@ window.MathJax = { - loader: {load: ['[tex]/ams']}, - tex: { - inlineMath: [["\\(", "\\)"]], - displayMath: [["\\[", "\\]"]], - processEscapes: true, - processEnvironments: true, - packages: {'[+]': ['ams']} - }, - options: { - ignoreHtmlClass: ".*|", - processHtmlClass: "arithmatex" - } - }; + loader: {load: ['[tex]/ams']}, + tex: { + inlineMath: [["\\(", "\\)"]], + displayMath: [["\\[", "\\]"]], + processEscapes: true, + processEnvironments: true, + packages: {'[+]': ['ams']} + }, + options: { + ignoreHtmlClass: ".*|", + processHtmlClass: "arithmatex" + } +}; + - document$.subscribe(() => { - MathJax.startup.output.clearCache() - MathJax.typesetClear() - MathJax.texReset() - MathJax.typesetPromise() - }) \ No newline at end of file +document$.subscribe(() => { + MathJax.startup.output.clearCache() + MathJax.typesetClear() + MathJax.texReset() + MathJax.typesetPromise() +}) From 93e37b2940800d2f5f5db5ec0eaa179477654075 Mon Sep 17 00:00:00 2001 From: raar1 Date: Wed, 4 Feb 2026 14:17:35 +0100 Subject: [PATCH 05/10] Add line from the mkdocs-jupyter repo config that seems to be necessary to make latex equations render in markdown cells --- mkdocs.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/mkdocs.yml b/mkdocs.yml index 3e91060..9221d03 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -73,6 +73,7 @@ plugins: - mkdocs-jupyter: execute: true include_source: True + custom_mathjax_url: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/latest.js?config=TeX-AMS_CHTML-full,Safe" # Styled blocks: https://squidfunk.github.io/mkdocs-material/reference/admonitions/#supported-types markdown_extensions: From 62589b2ac149a885b01cb00de79c3e3a68f5d53c Mon Sep 17 00:00:00 2001 From: raar1 Date: Wed, 4 Feb 2026 14:24:00 +0100 Subject: [PATCH 06/10] For now we will keep the examples/ dir in the root, and add a docs/examples dir with two notebooks to illustrate the rendering. T. This can be sorted out after Jorrit finishes the documentation PR --- .../additional_term.ipynb | 0 .../example_equation_formulation.ipynb | 0 .../example_exact_equation.ipynb | 0 ...xample_mibitrans_anatrans_comparison.ipynb | 387 +++++++ examples/example_walkthrough.ipynb | 967 ++++++++++++++++++ 5 files changed, 1354 insertions(+) rename {docs/examples => examples}/additional_term.ipynb (100%) rename {docs/examples => examples}/example_equation_formulation.ipynb (100%) rename {docs/examples => examples}/example_exact_equation.ipynb (100%) create mode 100644 examples/example_mibitrans_anatrans_comparison.ipynb create mode 100644 examples/example_walkthrough.ipynb diff --git a/docs/examples/additional_term.ipynb b/examples/additional_term.ipynb similarity index 100% rename from docs/examples/additional_term.ipynb rename to examples/additional_term.ipynb diff --git a/docs/examples/example_equation_formulation.ipynb b/examples/example_equation_formulation.ipynb similarity index 100% rename from docs/examples/example_equation_formulation.ipynb rename to examples/example_equation_formulation.ipynb diff --git a/docs/examples/example_exact_equation.ipynb b/examples/example_exact_equation.ipynb similarity index 100% rename from docs/examples/example_exact_equation.ipynb rename to examples/example_exact_equation.ipynb diff --git a/examples/example_mibitrans_anatrans_comparison.ipynb b/examples/example_mibitrans_anatrans_comparison.ipynb new file mode 100644 index 0000000..c4110de --- /dev/null +++ b/examples/example_mibitrans_anatrans_comparison.ipynb @@ -0,0 +1,387 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "initial_id", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from mibitrans.data.parameters import AttenuationParameters\n", + "from mibitrans.data.parameters import HydrologicalParameters\n", + "from mibitrans.data.parameters import ModelParameters\n", + "from mibitrans.data.parameters import SourceParameters\n", + "from mibitrans.transport.models import Anatrans\n", + "from mibitrans.transport.models import Mibitrans\n", + "from mibitrans.visualize.animation import animate_1d" + ] + }, + { + "cell_type": "markdown", + "id": "d24df86e0eb6973d", + "metadata": {}, + "source": [ + "### Exact solution versus untruncated solution\n", + "\n", + "\\begin{equation}\n", + "C(x,y,z,t) = \\sum_{i=1}^{n}\\left(C^*_{0,i}\\frac{x}{8\\sqrt{\\pi D^{'}_{x}}}\\exp(-\\gamma t)\n", + "\\cdot \\int_{0}^{t}\\left[\\frac{1}{\\tau^{\\frac{3}{2}}} \\exp\\left((\\gamma - \\lambda_{EFF})\\tau - \\frac{(x-v^{'}\\tau)^2}{4D^{'}_{x}\\tau}\\right)\n", + "\\cdot \\left\\{ERFC\\left(\\frac{y-Y_i}{2 \\sqrt{D^{'}_{y}\\tau}}\\right)-ERFC\\left(\\frac{y+Y_i}{2 \\sqrt{D^{'}_{y}\\tau}}\\right) \\right\\}\n", + "\\cdot \\left\\{ERFC\\left(\\frac{z-Z}{2 \\sqrt{D^{'}_{z}\\tau}}\\right)-ERFC\\left(\\frac{z+Z}{2 \\sqrt{D^{'}_{z}\\tau}}\\right) \\right\\}\\right] d\\tau \\right)\n", + "\\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29cec10a2865e2ca", + "metadata": {}, + "outputs": [], + "source": [ + "hydro = HydrologicalParameters(\n", + " velocity=0.277, # Flow velocity [m/d]\n", + " porosity=0.25, # Effective soil porosity [-]\n", + " alpha_x=10, # Longitudinal dispersivity, in [m]\n", + " alpha_y=1, # Transverse horizontal dispersivity, in [m]\n", + " alpha_z=0.1 # Transverse vertical dispersivity, in [m]\n", + ")\n", + "\n", + "att = AttenuationParameters(\n", + " retardation=1,\n", + " # Molecular diffusion, in [m2/day]\n", + " diffusion=0,\n", + " # Contaminant half life, in [days]\n", + " half_life=0#5*365\n", + ")\n", + "\n", + "source = SourceParameters(\n", + " source_zone_boundary=np.array([10]),\n", + " source_zone_concentration=np.array([11]),\n", + " depth=2.5,\n", + " total_mass=\"inf\"\n", + ")\n", + "\n", + "model = ModelParameters(\n", + " # Model extent in the longitudinal (x) direction in [m].\n", + " model_length = 800,\n", + " # Model extent in the transverse horizontal (y) direction in [m].\n", + " model_width = 150,\n", + " # Model duration in [days].\n", + " model_time = 5*365,\n", + " # Model grid discretization step size in the longitudinal (x) direction, in [m].\n", + " dx = 2,\n", + " # Model grid discretization step size in the transverse horizontal (y) direction, in [m].\n", + " dy = 1,\n", + " # Model time discretization step size, in [days]\n", + " dt = 25\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "54e9e7e491c986fe", + "metadata": {}, + "outputs": [], + "source": [ + "mbt_object = Mibitrans(hydro, att, source, model)\n", + "mbt_object.run()\n", + "\n", + "ana_object = Anatrans(hydro, att, source, model)\n", + "ana_object.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e323802bcb322adf", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3d53214555ede3a", + "metadata": {}, + "outputs": [], + "source": [ + "ani = animate_1d(x_axis_parameter=mbt_object.y,\n", + " y_axis_parameter=[mbt_object.cxyt[:,:,25], mbt_object.cxyt[:,:,50],\n", + " mbt_object.cxyt[:,:,100],mbt_object.cxyt[:,:,200],\n", + " ana_object.cxyt[:,:,25], ana_object.cxyt[:,:,50],\n", + " ana_object.cxyt[:,:,100], ana_object.cxyt[:,:,200]],\n", + " time_parameter=mbt_object.t,\n", + " y_colors=[\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\",\n", + " \"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\"],\n", + " y_names=[\"Mibitrans x=50m\", \"Mibitrans x=100m\", \"Mibitrans x=200m\", \"Mibitrans x=400m\",\n", + " \"Anatrans x=50m\", \"Anatrans x=100m\", \"Anatrans x=200m\", \"Anatrans x=400m\"],\n", + " linestyle=[\"-\", \"-\", \"-\", \"-\", \"--\", \"--\", \"--\", \"--\"])\n", + "plt.xlabel(\"y-position [m]\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78a4f70e30d859ea", + "metadata": {}, + "outputs": [], + "source": [ + "ani1 = animate_1d(x_axis_parameter=mbt_object.x,\n", + " y_axis_parameter=[mbt_object.cxyt[:,75,:], mbt_object.cxyt[:,84,:],\n", + " mbt_object.cxyt[:,100,:], mbt_object.cxyt[:,125,:],\n", + " ana_object.cxyt[:,75,:], ana_object.cxyt[:,84,:],\n", + " ana_object.cxyt[:,100,:], ana_object.cxyt[:,125,:]],\n", + " time_parameter=mbt_object.t,\n", + " y_colors=[\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\",\n", + " \"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\"],\n", + " y_names=[\"Mibitrans y=0m\", \"Mibitrans y=9m\", \"Mibitrans y=25m\", \"Mibitrans y=50m\",\n", + " \"Anatrans y=0m\", \"Anatrans y=9m\", \"Anatrans y=25m\", \"Anatrans y=50m\"],\n", + " linestyle=[\"-\", \"-\", \"-\", \"-\", \"--\", \"--\", \"--\", \"--\"])\n", + "plt.xlabel(\"x-position [m]\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "77b3c43e8ed23283", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "mbt_object.breakthrough(50, 0, color=\"darkgreen\", linestyle=\"-\", label=\"Mibitrans 50m\")\n", + "mbt_object.breakthrough(100, 0, color=\"limegreen\", linestyle=\"-\", label=\"Mibitrans 100m\")\n", + "mbt_object.breakthrough(200, 0, color=\"cornflowerblue\", linestyle=\"-\", label=\"Mibitrans 200m\")\n", + "mbt_object.breakthrough(400, 0, color=\"blue\", linestyle=\"-\", label=\"Mibitrans 400m\")\n", + "ana_object.breakthrough(50, 0, color=\"green\", linestyle=\"--\", label=\"Anatrans 50m\")\n", + "ana_object.breakthrough(100, 0, color=\"greenyellow\", linestyle=\"--\", label=\"Anatrans 100m\")\n", + "ana_object.breakthrough(200, 0, color=\"darkturquoise\", linestyle=\"--\", label=\"Anatrans 200m\")\n", + "ana_object.breakthrough(400, 0, color=\"dodgerblue\", linestyle=\"--\", label=\"Anatrans 400m\")\n", + "plt.title(\"Breakthrough plot of Mibitrans and Anatrans model, at various x locations.\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c38f7396fbd67f73", + "metadata": {}, + "source": [ + "#### Comparing model differences for various flow velocities and dispersivities" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae416dc23dd6b334", + "metadata": {}, + "outputs": [], + "source": [ + "parameter_list = [0.1, 0.5, 1, 5, 10]\n", + "output_list_mbt = []\n", + "output_list_ana = []\n", + "for par in parameter_list:\n", + " print(f\"starting run with par = {par}\")\n", + " mbt_object = Mibitrans(hydro, att, source, model)\n", + " mbt_object.hydrological_parameters.alpha_x = par\n", + " mbt_object.hydrological_parameters.alpha_y = par / 10\n", + " mbt_object.hydrological_parameters.alpha_z = par / 100\n", + " mbt_object.run()\n", + " output_list_mbt.append(mbt_object.cxyt)\n", + "\n", + " ana_object = Anatrans(hydro, att, source, model)\n", + " ana_object.hydrological_parameters.alpha_x = par\n", + " ana_object.hydrological_parameters.alpha_y = par / 10\n", + " ana_object.hydrological_parameters.alpha_z = par / 100\n", + " ana_object.run()\n", + " output_list_ana.append(ana_object.cxyt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93a0274d7116f5f", + "metadata": {}, + "outputs": [], + "source": [ + "colors_mbt = [\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\", \"indigo\"]\n", + "colors_ana = [\"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\", \"darkviolet\"]\n", + "plt.figure()\n", + "for i, conc in enumerate(output_list_mbt):\n", + " plt.plot(mbt_object.x, conc[-1,75,:], color=colors_mbt[i], label=f\"Mibitrans a={parameter_list[i]}m\")\n", + "for i, conc in enumerate(output_list_ana):\n", + " plt.plot(ana_object.x, conc[-1,75,:], color=colors_ana[i], linestyle=\"--\", label=f\"Anatrans a={parameter_list[i]}m\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea1d7ea112addeff", + "metadata": {}, + "outputs": [], + "source": [ + "colors_mbt = [\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\", \"indigo\"]\n", + "colors_ana = [\"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\", \"darkviolet\"]\n", + "plt.figure()\n", + "for i, conc in enumerate(output_list_mbt):\n", + " plt.plot(mbt_object.y, conc[-1,:,175], color=colors_mbt[i],\n", + " label=f\"Mibitrans a={parameter_list[i]}m\")\n", + "for i, conc in enumerate(output_list_ana):\n", + " plt.plot(ana_object.y, conc[-1,:,175], color=colors_ana[i],\n", + " linestyle=\"--\", label=f\"Anatrans a={parameter_list[i]}m\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e8640aff2864ca7", + "metadata": {}, + "outputs": [], + "source": [ + "parameter_list_v = [0.1, 0.1, 0.1, 0.1, 0.1]\n", + "parameter_list_a = [10, 5, 2, 1, 0.5]\n", + "output_list_mbt_v = []\n", + "output_list_ana_v = []\n", + "for i in range(len(parameter_list_v)):\n", + " parv = parameter_list_v[i]\n", + " para = parameter_list_a[i]\n", + " print(f\"starting run with parv = {parv} and para = {para}\")\n", + " mbt_object = Mibitrans(hydro, att, source, model)\n", + " mbt_object.hydrological_parameters.velocity = parv\n", + " mbt_object.hydrological_parameters.alpha_x = para\n", + " mbt_object.hydrological_parameters.alpha_y = para / 10\n", + " mbt_object.hydrological_parameters.alpha_z = 0#para / 100\n", + " mbt_object.run()\n", + " output_list_mbt_v.append(mbt_object.cxyt)\n", + "\n", + " ana_object = Anatrans(hydro, att, source, model)\n", + " ana_object.hydrological_parameters.velocity = parv\n", + " ana_object.hydrological_parameters.alpha_x = para\n", + " ana_object.hydrological_parameters.alpha_y = para / 10\n", + " ana_object.hydrological_parameters.alpha_z = 0#para / 100\n", + " ana_object.run()\n", + " output_list_ana_v.append(ana_object.cxyt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5383ae968db9bf6", + "metadata": {}, + "outputs": [], + "source": [ + "colors_mbt = [\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\", \"indigo\"]\n", + "colors_ana = [\"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\", \"darkviolet\"]\n", + "plt.figure()\n", + "for i, conc in enumerate(output_list_mbt_v):\n", + " plt.plot(mbt_object.x, conc[-1,75,:], color=colors_mbt[i],\n", + " label=f\"Mbt v={parameter_list_v[i]}m/d, a={parameter_list_a[i]}m\")\n", + "for i, conc in enumerate(output_list_ana_v):\n", + " plt.plot(ana_object.x, conc[-1,75,:], color=colors_ana[i],\n", + " linestyle=\"--\", label=f\"Ana v={parameter_list_v[i]}m/d, a={parameter_list_a[i]}m\")\n", + "plt.xlim((0,300))\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f38d34e6b66ce99e", + "metadata": {}, + "outputs": [], + "source": [ + "colors_mbt = [\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\", \"indigo\"]\n", + "colors_ana = [\"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\", \"darkviolet\"]\n", + "plt.figure()\n", + "for i, conc in enumerate(output_list_mbt_v):\n", + " plt.plot(mbt_object.y, conc[-1, :, 75], color=colors_mbt[i], label=f\"Mibitrans a={parameter_list_v[i]}m\")\n", + "for i, conc in enumerate(output_list_ana_v):\n", + " plt.plot(ana_object.y, conc[-1, :, 75], color=colors_ana[i], linestyle=\"--\",\n", + " label=f\"Anatrans a={parameter_list[i]}m\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f01c7d5012b64963", + "metadata": {}, + "outputs": [], + "source": [ + "mbtv=output_list_mbt_v\n", + "anav=output_list_ana_v\n", + "ani = animate_1d(x_axis_parameter=mbt_object.x,\n", + " y_axis_parameter=[mbtv[0][:,75,:], mbtv[1][:,75,:],mbtv[2][:,75,:],\n", + " mbtv[3][:,75,:], mbtv[4][:,75,:],\n", + " anav[0][:,75,:], anav[1][:,75,:], anav[2][:,75,:],\n", + " anav[3][:,75,:], anav[4][:,75,:]],\n", + " time_parameter=mbt_object.t,\n", + " y_colors=[\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\", \"indigo\",\n", + " \"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\", \"violet\"],\n", + " y_names=[\"Mibitrans a=0.1\", \"Mibitrans a=0.5\", \"Mibitrans a=1\", \"Mibitrans a=5\", \"Mibitrans a=10\",\n", + " \"Anatrans a=0.1\", \"Anatrans a=0.5\", \"Anatrans a=1\", \"Anatrans a=5\", \"Anatrans a=10\"],\n", + " linestyle=[\"-\", \"-\", \"-\", \"-\", \"-\", \"--\", \"--\", \"--\", \"--\", \"--\"])\n", + "plt.xlabel(\"y-position [m]\")\n", + "plt.xlim((0,350))\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f70026152c8958e", + "metadata": {}, + "outputs": [], + "source": [ + "mbtv=output_list_mbt_v\n", + "anav=output_list_ana_v\n", + "ani = animate_1d(x_axis_parameter=mbt_object.y,\n", + " y_axis_parameter=[mbtv[0][:,:,50], mbtv[1][:,:,50], mbtv[2][:,:,50],\n", + " mbtv[3][:,:,50], mbtv[4][:,:,50],\n", + " anav[0][:,:,50], anav[1][:,:,50], anav[2][:,:,50],\n", + " anav[3][:,:,50], anav[4][:,:,50]],\n", + " time_parameter=mbt_object.t,\n", + " y_colors=[\"darkgreen\", \"limegreen\", \"cornflowerblue\", \"blue\", \"indigo\",\n", + " \"green\", \"greenyellow\", \"darkturquoise\", \"dodgerblue\", \"violet\"],\n", + " y_names=[\"Mibitrans a=0.1\", \"Mibitrans a=0.5\", \"Mibitrans a=1\", \"Mibitrans a=5\", \"Mibitrans a=10\",\n", + " \"Anatrans a=0.1\", \"Anatrans a=0.5\", \"Anatrans a=1\", \"Anatrans a=5\", \"Anatrans a=10\"],\n", + " linestyle=[\"-\", \"-\", \"-\", \"-\", \"-\", \"--\", \"--\", \"--\", \"--\", \"--\"])\n", + "plt.xlabel(\"y-position [m]\")\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/example_walkthrough.ipynb b/examples/example_walkthrough.ipynb new file mode 100644 index 0000000..940cbaf --- /dev/null +++ b/examples/example_walkthrough.ipynb @@ -0,0 +1,967 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from mibitrans.data.parameter_information import ElectronAcceptors\n", + "from mibitrans.data.parameter_information import UtilizationFactor\n", + "from mibitrans.data.parameters import AttenuationParameters\n", + "from mibitrans.data.parameters import HydrologicalParameters\n", + "from mibitrans.data.parameters import ModelParameters\n", + "from mibitrans.data.parameters import SourceParameters\n", + "from mibitrans.transport.models import Anatrans\n", + "from mibitrans.transport.models import Bioscreen\n", + "from mibitrans.transport.models import Mibitrans\n", + "from mibitrans.visualize.plot_line import breakthrough\n", + "from mibitrans.visualize.plot_line import centerline\n", + "from mibitrans.visualize.plot_line import transverse\n", + "from mibitrans.visualize.plot_surface import plume_2d\n", + "from mibitrans.visualize.plot_surface import plume_3d" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The showcase of mibitrans below uses example data originating from BIOSCREEN, the Excel based modeling software this package took inspiration from. Describes field conditions on the Keesler Air Force Base in Mississippi, USA. Movement of mixed BTEX plume from 1989 to 1995. For more details, see Newell et al. (1997)\n", + "\n", + "As it was developed in the USA, length units are in ft. We'll use a conversion factor to express length in m instead.\n", + "\n", + "Newell, C. J., McLeod, R. K., & Gonzales, J. R. (1997). BIOSCREEN natural attenuation decision support\n", + "system version 1.4 revisions, Tech. rep., U.S. EPA." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ft = 3.281 # factor to convert ft to m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Input by dataclasses\n", + "\n", + "mibitrans use dataclasses located in mibitrans.data.read to handle data input. This can be There are five input dataclasses, each for a different category of parameters. To avoid any mistakes, units of input parameters should be the ones specified. When using different units, make sure that they are consistent throughout the entire modelling process." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Hydrological parameters\n", + "\n", + "Contains parameters that are inherent to the aquifer properties; flow velocity, porosity and dispersivity. Flow velocity can alternatively be calculated from the hydraulic conductivity and hydraulic gradient. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "hydro = HydrologicalParameters(\n", + " velocity=113.8/ft/365, # Groundwater flow velocity, in [m/day]\n", + " porosity=0.25, # Effective soil porosity [-]\n", + " alpha_x=13.3/ft, # Longitudinal dispersivity, in [m]\n", + " alpha_y=1.3/ft, # Transverse horizontal dispersivity, in [m]\n", + " alpha_z=0 # Transverse vertical dispersivity, in [m]\n", + ")\n", + "\n", + "# Alternative by specifying hydraulic gradient and hydraulic conductivity\n", + "hydro = HydrologicalParameters(\n", + " h_gradient=0.048, # Hydraulic gradient [-]\n", + " h_conductivity=0.495, # Hydraulic conductivity [m/day]\n", + " porosity=0.25, # Effective soil porosity [-]\n", + " alpha_x=13.3/ft, # Longitudinal dispersivity, in [m]\n", + " alpha_y=1.3/ft, # Transverse horizontal dispersivity, in [m]\n", + " alpha_z=0 # Transverse vertical dispersivity, in [m]\n", + ")\n", + "\n", + "# And then check the value:\n", + "print(\"The calculated groundwater flow velocity is: \", hydro.velocity, r\"m/d\")\n", + "\n", + "# Any other input parameters can be requested by specifying the argument;\n", + "print(f\"The dispersivity values are: {hydro.alpha_x}m, {hydro.alpha_y}m and {hydro.alpha_z}m,\"\n", + " \" for the x, y and z directions respectively.\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Attenuation parameters\n", + "Handles all parameters related to adsorption, diffusion and degradation.\n", + "\n", + "##### Attenuation: retardation and diffusion\n", + "Retardation can be simply given as a value >= 1. Alternatively, the adsorption is calculated from the soil bulk density, paratition coefficient and the fraction of organic carbon in the soil. Note that calculation of retardation factor requires porosity as well, which is already provided in HydrologicalParameters. The calculation of retardation will therefore be automatically performed in the analytical equation. It can be manually calculated using the calculate_retardation method as well.\n", + "\n", + "Here, molecular diffusion can be given as well. It is 0 by default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "att = AttenuationParameters(\n", + " # Retardation factor for transported contaminants\n", + " retardation=1\n", + ")\n", + "\n", + "# Alternatively, calculate the retardation factor by supplying soil and contaminant properties\n", + "att = AttenuationParameters(\n", + " # Soil bulk density in [g/m^3]\n", + " bulk_density=1.7,\n", + " # Partition coefficient of the transported contaminant to soil organic matter, in [m^3/g]\n", + " partition_coefficient=38,\n", + " # Fraction of organic material in the soil [-]\n", + " fraction_organic_carbon=5.7e-5\n", + ")\n", + "\n", + "# Calculate the retardation factor beforehand to see its value by specifying the porosity\n", + "att.calculate_retardation(porosity=hydro.porosity)\n", + "\n", + "# And then check the value:\n", + "print(\"The calculated retardation value is: \", att.retardation)\n", + "\n", + "att = AttenuationParameters(\n", + " # Retardation factor for transported contaminants\n", + " retardation=1,\n", + " # Molecular diffusion, in [m2/day]\n", + " diffusion=1e-5\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Attenuation: degradation parameters\n", + "\n", + "Linear decay models only need either the contaminant decay rate or half life. Input for the instant reaction model is somewhat more involved, and is therefore done with a class method of the model classes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "att = AttenuationParameters(\n", + " # Contaminant first order decay rate in [1/days]\n", + " decay_rate=0.0127\n", + ")\n", + "\n", + "# Alternatively, specify the contaminant half life\n", + "att = AttenuationParameters(\n", + " # Contaminant half life, in [days]\n", + " half_life=54.75,\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the examples below, we'll set the attenuation parameters to the desired settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "att = AttenuationParameters(\n", + " # Soil bulk density in [g/m^3]\n", + " bulk_density=1.7,\n", + " # Partition coefficient of the transported contaminant to soil organic matter, in [m^3/g]\n", + " partition_coefficient=38,\n", + " # Fraction of organic material in the soil [-]\n", + " fraction_organic_carbon=5.7e-5,\n", + " # Molecular diffusion, in [m2/day]\n", + " diffusion=0,\n", + " # Contaminant half life, in [days]\n", + " half_life=365#54.75,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Source parameters\n", + "\n", + "Takes input of the dimensions of and concentrations at the contaminant source. The source is treated as a seperate phase, which dissolves into the groundwater over time. The source is assumed to be in symmetrical in its center, and concentrations decrease from the center to the fringes. Furthermore, the source is assumed to be constant over its depth. The transverse horizontal dimension (or width) of the source is divided into zones, which span a certain distance measured from the center of the source, each with an associated concentration. For example, a source can have a concentration of $10g/m^3$ $7m$ left and right from the source center, and a concentration of $5g/m^3$ up to $20m$ from the source center. This would then be entered as source_zone_boundary = $[7,20]$ and source_zone_concentration = $[10,5]$. The source can be a single zone with a single concentration as well.\n", + "\n", + "By giving a total mass of the contaminant source, the amount of solid-phase contaminant, and with that, the source zone concentrations, diminish over time. The rate at which this occurs depends on the flow velocity in the aquifer and the size of the source zone. The source zone can be set to be considered infinite as well, meaning that the concentrations will not diminish over time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Input for a simple source zone with a width of 19.8m (65ft) and a continuous input (infinite source mass)\n", + "source = SourceParameters(\n", + " # Source zone boundaries, in [m] (simply using a float instead of a numpy array for\n", + " # single source zone input will work as well)\n", + " source_zone_boundary=np.array([65/ft]),\n", + " # Source zone concentrations, in [g/m^3]\n", + " source_zone_concentration=np.array([5]),\n", + " # Source depth extent, in [m]\n", + " depth=10/ft,\n", + " # Source mass, considered infinite\n", + " total_mass=\"inf\"\n", + ")\n", + "\n", + "# Alternatively, specify a source mass to allow for source decay\n", + "source = SourceParameters(\n", + " source_zone_boundary=np.array([7/ft, 37/ft, 65/ft]),\n", + " source_zone_concentration=np.array([13.68, 2.508, 0.057]),\n", + " depth=10/ft,\n", + " total_mass=2000000\n", + ")\n", + "\n", + "# Visualize what the source zone looks like to check your input:\n", + "source.visualize()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Model parameters\n", + "\n", + "Accepts input for the model dimensions and discretization. Model length is the extent in the (x) direction parallel to the groundwater flow direction. The model width is the extent of the model perpendicular (y) to the groundwater flow direction. Step size of the spatial dimensions is handled with dx and dy. Ensure that the source zone fits inside of the given model width. If step sizes are not given, a ratio of model_length (1/100), model_width (1/50) and model_time (1/10) is used by default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model = ModelParameters(\n", + " # Model extent in the longitudinal (x) direction in [m].\n", + " model_length = 320/ft,\n", + " # Model extent in the transverse horizontal (y) direction in [m].\n", + " model_width = 100/ft,\n", + " # Model duration in [days].\n", + " model_time = 6 * 365,\n", + " # Model grid discretization step size in the longitudinal (x) direction, in [m].\n", + " dx = 1/ft,\n", + " # Model grid discretization step size in the transverse horizontal (y) direction, in [m].\n", + " dy = 1/ft,\n", + " # Model time discretization step size, in [days]\n", + " dt = 365 / 5\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Input checking\n", + "\n", + "The dataclass inputs evaluates if the required input parameters are present, if they are of the correct data type and if they are in the expected domain." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# If not all required parameters are specified, an error will be shown;\n", + "# HydrologicalParameters needs the porosity, longitudinal dispersivity and transverse horizontal dispersivity as well.\n", + "fake_hydro = HydrologicalParameters(velocity = 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# If the input datatype is not what is expected; source_zone_concentration expects a numpy array\n", + "# of the same length as the array given in source_zone_boundary.\n", + "fake_source = SourceParameters(\n", + " source_zone_boundary = np.array([1,2,3]),\n", + " source_zone_concentration = \"this is a string, not an array\",\n", + " depth = 10,\n", + " total_mass = \"inf\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Same goes for if the input parameter has a value outside its valid domain; retardation should have a value >= 1\n", + "fake_att = AttenuationParameters(\n", + " retardation = 0.1\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Analytical models\n", + "\n", + "The various models located in mibitrans.transport.models. Currently, three models are implemented; 'Mibitrans', 'Anatrans' and 'Bioscreen'. Each with a distinct analytical solution, which are introduced below.\n", + "\n", + "Input to the model is provided through the dataclasses showcased above. Each input argument has the same name as the dataclass, but in snake_case instead of PascalCase (e.g. hydrological_paramaters takes in the data class HydrologicalParameters), to follow Python convention. The models are implemented as classes, which are first initialized through making the object. This will generate the model grid and do some prior parameter calculations. The model object can then be used to call attributes, or class methods. After initialization, no results have been calculated yet. For this, use the 'run' class method. This generates a results class object, which contains the concentration distribution as a 3D numpy array (indexed as $C[t,y,x]$) and the parameters used for this specific run." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Mibitrans model\n", + "\n", + "The Mibitrans model class uses the exact analytical solution implemented by Karanovic (2007) in the Excel based BIOSCREEN-AT, and added source depletion, akin to that implemented in its predecessor BIOSCREEN by Newell et al. (1997). This model is based on the Wexler (1992) solution. The Mibitrans model allows for the same method as used in BIOSCREEN-AT, but expands it by allowing multiple source zones (by means of superposition) and including the instant reaction model. These were present in the original BIOSCREEN, but not reimplemented in BIOSCREEN-AT. Using a single source zone in this model, and not using the instant reaction option will make the Mibitrans solution resolve to the equation described in Karanovic (2007). Which in turn resolves to the Wexler (1992) solution if source depletion is disabled.\n", + "\n", + "As the namesake model of this package, it is the recommended model to use. The other models introduce a margin of error by making some assumptions. However, since this model requires evaluation of an integral, computation time might be longer, depending on model discretization. The exact nature of these differences is too much to go into detail here, but is elaborated upon in the theoretical background.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Initializing the model object\n", + "mbt_object = Mibitrans(\n", + " hydrological_parameters=hydro,\n", + " attenuation_parameters=att,\n", + " source_parameters=source,\n", + " model_parameters=model\n", + ")\n", + "\n", + "# Can use model object to request various attributes\n", + "print(\"x discretization: \", mbt_object.x)\n", + "print(\"Retarded flow velocity: \", mbt_object.rv)\n", + "print(\"x steps\", len(mbt_object.x), \"y steps\", len(mbt_object.y), \"t steps\", len(mbt_object.t))\n", + "\n", + "# Run the model once initialized and obtain the results\n", + "mbt_results = mbt_object.run()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "?Mibitrans" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Once the model has run, results are contained in the cxyt attribute\n", + "model_cxyt = mbt_object.cxyt\n", + "\n", + "# cxyt is indexed as [time, y-position, x-position]\n", + "# Thus to get the concentration at the last time step, in the center of the plume for all x:\n", + "plume_center = model_cxyt[-1, 132//2, :]\n", + "\n", + "# mibitrans has build-in visualization methods (see visualization section for more details).\n", + "mbt_results.centerline()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To perform the model with biodegradation modelled as an instant reaction, use the instant_reaction class method. Input for the instant reaction model is somewhat more involved, needing electron donor and acceptor concentrations. For utilization factors (amount of electron donor/acceptor used/generated by biodegradation), the values for BTEX degradation are used by default, but custom values can be given. For more specifics about the underlying principles and assumptions, see the theory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# For streamlined input, use the ElectronAcceptor dataclass;\n", + "ea = ElectronAcceptors(\n", + " # Difference between background oxygen and current oxygen concentration in groundwater, in [g/m^3]\n", + " delta_oxygen=1.65,\n", + " # Difference between background nitrate and current nitrate concentration in groundwater, in [g/m^3]\n", + " delta_nitrate=0.7,\n", + " # Current ferrous iron concentration in groundwater, in [g/m^3]\n", + " ferrous_iron=16.6,\n", + " # Difference between background sulfate and current sulfate concentration in groundwater, in [g/m^3]\n", + " delta_sulfate=22.4,\n", + " # Current methane concentration in groundwater, in [g/m^3]\n", + " methane=6.6\n", + ")\n", + "\n", + "mbt_object.instant_reaction(electron_acceptors=ea)\n", + "\n", + "# Can adapt utilization factors if needed\n", + "uf = UtilizationFactor(\n", + " # utilization factor of oxygen, as mass of oxygen consumed per mass of biodegraded contaminant [g/g].\n", + " util_oxygen=2,\n", + " # utilization factor of nitrate, as mass of nitrate consumed per mass of biodegraded contaminant [g/g].\n", + " util_nitrate=1,\n", + " # utilization factor of ferrous iron, as mass of ferrous iron generated per mass of biodegraded contaminant [g/g].\n", + " util_ferrous_iron=4,\n", + " # utilization factor of sulfate, as mass of sulfate consumed per mass of biodegraded contaminant [g/g].\n", + " util_sulfate=3,\n", + " # utilization factor of methane, as mass of methane generated per mass of biodegraded contaminant [g/g].\n", + " util_methane=5,\n", + ")\n", + "\n", + "mbt_object.instant_reaction(electron_acceptors=ea, utilization_factor=uf)\n", + "\n", + "# Alternatively, electron acceptors and/or utilization factors can be entered as a dictionary\n", + "mbt_object.instant_reaction(\n", + " electron_acceptors={\n", + " # Difference between background oxygen and current oxygen concentration in groundwater, in [g/m^3]\n", + " \"delta_oxygen\":1.65,\n", + " # Difference between background nitrate and current nitrate concentration in groundwater, in [g/m^3]\n", + " \"delta_nitrate\":0.7,\n", + " # Current ferrous iron concentration in groundwater, in [g/m^3]\n", + " \"ferrous_iron\":16.6,\n", + " # Difference between background sulfate and current sulfate concentration in groundwater, in [g/m^3]\n", + " \"delta_sulfate\":22.4,\n", + " # Current methane concentration in groundwater, in [g/m^3]\n", + " \"methane\":6.6,\n", + " }\n", + ")\n", + "print(mbt_object.biodegradation_capacity)\n", + "# Note that using instant_reaction also resets the utilization factor to default.\n", + "print(\"electron acceptor concentrations: \", mbt_object.electron_acceptors)\n", + "print(\"electron acceptor utilization factors: \", mbt_object.utilization_factor)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Be aware that electron acceptor concentrations and utilization factors can only be changed using the instant_reaction method. Changing the properties directly will not work.\n", + "\n", + "Now that instant reaction parameters are provided, the model can be run:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mbt_instant_results = mbt_object.run()\n", + "\n", + "mbt_instant_results.centerline()\n", + "plt.show()\n", + "\n", + "# Switch between instant reaction and linear model\n", + "mbt_object.mode = \"linear\"\n", + "print(\"The current model mode is: \", mbt_object.mode)\n", + "# Running now will show the results of the linear model\n", + "# Switch back to instant reaction as follows:\n", + "mbt_object.mode = \"instant_reaction\"\n", + "print(\"The current model mode is: \", mbt_object.mode)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use the sample method to get the concentration at a specific location and time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Instead, sample a specific location at any point along the plume and any point in time\n", + "concentration = mbt_object.sample(\n", + " x_position=100,\n", + " y_position=0,\n", + " time=10*365)\n", + "print(\"The concentration at 150m downstream from the source after 10 years is:\", concentration, r\"g/m^3\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Using the verbose flag, integration steps are printed to console.\n", + "# Usefull for longer runs to track progress\n", + "mbt_object = Mibitrans(\n", + " hydrological_parameters=hydro,\n", + " attenuation_parameters=att,\n", + " source_parameters=source,\n", + " model_parameters=model,\n", + " verbose=True\n", + ")\n", + "mbt_results = mbt_object.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Other models\n", + "\n", + "The two models other than the Mibitrans model share the same functionalities and properties as showcase above. The only difference is the implementation of calculation.\n", + "\n", + "#### Anatrans model\n", + "The equation used for the Anatrans model has the assumption that C(x,y,z,t) = C(x,t) * C(y,t) * C(z,t). Then, the 3D ADE can be broken up in three separate differential equations which can be solved individually. For C(x,t) the solution is given in Bear (1979), C(y,t) and C(z,t) can be derived from Crank (1975). The Anatrans model is the combination of these solutions, with addition of source depletion, source superposition and instant reaction model, described in Newell et al. (1997). The solution of Newell et al. (1997) is based of the Domenico (1987) solution, a truncated version of the equation described above, which introduces an error with a size dependent on the ratio of flow velocity and longitudinal dispersivity. Anatrans instead uses the fully untruncated version." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ana_object = Anatrans(\n", + " hydrological_parameters=hydro,\n", + " attenuation_parameters=att,\n", + " source_parameters=source,\n", + " model_parameters=model\n", + ")\n", + "\n", + "ana_results = ana_object.run()\n", + "ana_results.centerline()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Bioscreen model\n", + "\n", + "This model is an exact implementation of the transport equations implemented in the BIOSCREEN screening model of Newell et al. (1997), which is based on the Domenico (1987) analytical model. Using a truncated version of the equation used in the Anatrans model. This model is implemented as a method of comparison with the original BIOSCREEN software. And is included for legacy reasons, since it is the first model implemented in the mibitrans package, serving as a basis for the other models. However, caution should be taken when using this model, since a varying error is introduced by using the truncated analytical solution. The error is most prominent for shorter times and distances from the source, and depends on the ratio of flow velocity and longitudinal dispersivity. For modelling, the Anatrans (untruncated approximate solution) and Mibitrans (exact analytical solution) models are recommended instead." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bio_object = Bioscreen(\n", + " hydrological_parameters=hydro,\n", + " attenuation_parameters=att,\n", + " source_parameters=source,\n", + " model_parameters=model\n", + ")\n", + "bio_results = bio_object.run()\n", + "bio_results.centerline()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The model class check if they receive all required and the correct input dataclasses, ensuring model calculations will be performed without error." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Will not work; hydrological_parameters should be a HydrologicalParameters class object.\n", + "fake_mbt = Mibitrans(\n", + " hydrological_parameters=model,\n", + " attenuation_parameters=att,\n", + " source_parameters=source,\n", + " model_parameters=model\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization\n", + "\n", + "Mibitrans has various ways to visualize the model results. Above the centerline plot was already showcased, which plots the concentration distribution over the center of the contaminant plume. Plotting can be performed either as class method (i.e. model_object.centerline), or by calling the centerline function specifically; (i.e. centerline(model_object). The 1D plotting functions are located in mibitrans.visualize.plot_line. Multidimensional plots are located in mibitrans.visualize.plot_surface." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Lets make some different model objects to use in the plotting, takes a couple of seconds to run\n", + "\n", + "#Mibitrans model\n", + "mbt_object = Mibitrans(hydro, att, source, model)\n", + "mbt_lineardecay = mbt_object.run()\n", + "\n", + "mbt_object.attenuation_parameters.decay_rate = 0\n", + "mbt_nodecay = mbt_object.run()\n", + "\n", + "mbt_object.instant_reaction(electron_acceptors=ea)\n", + "mbt_instant = mbt_object.run()\n", + "\n", + "#Anatrans model\n", + "ana_object = Anatrans(hydro, att, source, model)\n", + "ana_lineardecay = ana_object.run()\n", + "\n", + "ana_object.attenuation_parameters.decay_rate = 0\n", + "ana_nodecay = ana_object.run()\n", + "\n", + "ana_object.instant_reaction(electron_acceptors=ea)\n", + "ana_instant = ana_object.run()\n", + "\n", + "#Bioscreen model\n", + "bio_object = Bioscreen(hydro, att, source, model)\n", + "bio_lineardecay = bio_object.run()\n", + "\n", + "bio_object.attenuation_parameters.decay_rate = 0\n", + "bio_nodecay = bio_object.run()\n", + "\n", + "bio_object.instant_reaction(electron_acceptors=ea)\n", + "bio_instant = bio_object.run()\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plotting centerline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pass model object to plotting function, by default, the last time step is used.\n", + "centerline(mbt_nodecay)\n", + "plt.show()\n", + "\n", + "# If you want to plot somewhere and sometime specific, use the time and y_position arguments.\n", + "# It gives the concentration profile at the step closest to what you specified.\n", + "centerline(mbt_nodecay, time=3*365, y_position=5)\n", + "#plt.title(\"No degradation Domenico model 5m away from center, t=6years\")\n", + "plt.show()\n", + "\n", + "# If you want you can change the plot settings to the ones you prefer\n", + "centerline(mbt_nodecay, time=6*365)\n", + "plt.title(\"Better title than the one generated automatically\")\n", + "plt.xlabel(\"I have changed!\")\n", + "plt.ylabel(\"And so have I\")\n", + "plt.show()\n", + "\n", + "legend = [\"no degradation\", \"linear decay\", \"instant reaction\"]\n", + "\n", + "# Instead of a single model, all line visualization functions accept a list of models to be displayed\n", + "# together in a single plot.\n", + "centerline([mbt_nodecay, mbt_lineardecay, mbt_instant], time=6*365, legend_names=legend)\n", + "plt.title(\"Mibitrans model at plume center for different models, t=6years\")\n", + "plt.show()\n", + "\n", + "# Keyword arguments for plt.plot can be passed on through the function\n", + "# Calling function separately per model gives more control over plot layout\n", + "centerline(mbt_nodecay, time=6*365, linestyle=\"--\", color=\"green\", label=legend[0])\n", + "centerline(mbt_lineardecay, time=6*365, linestyle=\"-.\", color=\"red\", label=legend[1])\n", + "plt.title(\"No and linear degradation Mibitrans model at plume center, t=6years\")\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "# Alternatively, use the inherent plotting modules of the model objects\n", + "mbt_nodecay.centerline(time = 365, linestyle=\"--\", color=\"green\", label=\"Mibitrans\")\n", + "ana_nodecay.centerline(time = 365, linestyle=\"-.\", color=\"blue\", label=\"Anatrans\")\n", + "bio_nodecay.centerline(time = 365, linestyle=\":\", color=\"red\", label=\"Bioscreen\")\n", + "plt.title(\"Concentration distribution for different contaminant transport models after 1 year.\")\n", + "plt.legend()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plotting transverse distribution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Concentration distribution can also be plotted in the transverse direction\n", + "transverse(ana_nodecay, x_position=80, time=6*365, linestyle=\"--\", color=\"green\", label=\"no degradation\")\n", + "transverse(ana_instant, x_position=80, time=6*365, linestyle=\"-.\", color=\"red\", label=\"instant reaction\")\n", + "\n", + "plt.title(\"No degradation and instant reaction Domenico model at x=80m, t=6years\")\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "# Or using the class methods\n", + "ana_lineardecay.transverse(x_position=80, time=6*365)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plotting breakthrough curve" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Concentration distribution can also be plotted in the transverse direction\n", + "breakthrough([bio_nodecay, bio_lineardecay, bio_instant], x_position=80,\n", + " legend_names=[\"no degradation\", \"linear decay\", \"instant reaction\"])\n", + "plt.title(\"Breakthrough curve of Bioscreen model for different degradation settings, at x=80m, t=6years\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plotting in 2D" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the x and y concentration distribution for the Mibitrans model, uses plt.pcolormesh\n", + "plume_2d(mbt_nodecay, time=6*365)\n", + "plt.title(\"Contaminant plume with no degradation Mibitrans model, t = 6 years\")\n", + "plt.show()\n", + "\n", + "# Function passes plt.colormesh keyword arguments\n", + "plume_2d(mbt_lineardecay, time=6*365, cmap=\"magma\")\n", + "plt.title(\"Contaminant plume with linear degradation Domenico model, t = 6 years\")\n", + "plt.show()\n", + "\n", + "# Once again also can be accessed through the class method\n", + "mbt_instant.plume_2d(time=4*365, cmap=\"plasma\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plotting 2D in 3D" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the x and y concentration distribution for no degradation decay model, uses plot_surface\n", + "plume_3d(ana_nodecay, time=6*365)\n", + "plt.title(\"Contaminant plume with no degradation Anatrans model, t = 6 years\")\n", + "plt.show()\n", + "\n", + "# Function passes plot_surface keyword arguments\n", + "plume_3d(ana_lineardecay, time=6*365, cmap=\"viridis\")\n", + "plt.title(\"Contaminant plume with linear degradation Anatrans model, t = 6 years\")\n", + "plt.show()\n", + "\n", + "# Function returns 'ax' object, use this to change view point of plot\n", + "# And can be accessed through class methods\n", + "ax = ana_instant.plume_3d(time=6*365, cmap=\"viridis\")\n", + "plt.title(\"Contaminant plume with linear degradation Anatrans model, t = 6 years\")\n", + "ax.view_init(elev=15, azim=340)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Animate plots\n", + "\n", + "All plots mentioned above in the visualization section have the option to be animated, which also can be saved as a file. Multiple models can be combined in a single animation. Make sure that parameters passed to the functions are inside the domain of all models. As each animation frame is a model time step, all models should have the exact same dt, otherwise, the animation will not show the correct temporal change in concentration." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Needed to show animations in Jupyter Notebooks\n", + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Output needs to be assigned to variable for animation to work\n", + "ani = centerline([mbt_nodecay, mbt_lineardecay, mbt_instant], time=6*365, legend_names=legend, animate=True)\n", + "plt.show()\n", + "\n", + "# Animation of breakthrough curve instead shows a timelapse of drawing of each curve\n", + "ani1 = breakthrough([bio_nodecay, bio_lineardecay, bio_instant], x_position=20, legend_names=legend, animate=True)\n", + "plt.show()\n", + "\n", + "# For the 3d surface plot, an entirely new plot needs to be generated per time step\n", + "# This can cause slightly longer execution times\n", + "ani2 = ana_instant.plume_3d(time=6*365, animate=True, cmap=\"viridis\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save animation as .gif (or other desired file format), this may take a while\n", + "# Adjust animation speed with optional fps argument\n", + "ani2.save(\"plume_animation.gif\", fps=10)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mass balance\n", + "\n", + "To gain numerical information about mass transport in the model area, use the mass balance method of the Results object. Each of the mass balance elements is a class property. Descriptions of each mass balance element can be found in analysis.mass_balance.MassBalance." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up a mass balance object of the Mibitrans instant reaction model.\n", + "mb = mbt_instant.mass_balance()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As the warning shows, part of the plume extents beyond the model boundary, which makes the mass balance inaccurate. Increasing model length will remediate this." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mbt_object.model_parameters.model_length = 500 / ft\n", + "results = mbt_object.run()\n", + "mb = results.mass_balance()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The entire plume is now inside the model extent." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results.centerline()\n", + "plt.show()\n", + "results.transverse(20, label=\"x = 20m\")\n", + "results.transverse(60, label=\"x = 60m\")\n", + "results.transverse(100, label=\"x = 100m\")\n", + "results.transverse(140, label=\"x = 140m\")\n", + "plt.legend()\n", + "plt.title(\"Transverse plot of Mibitrans instant reaction model, at t=6 years\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# By default, mass balance gives results for each model time step\n", + "print(mb.plume_mass)\n", + "# And then can easily be plotted over time\n", + "plt.plot(mb.t, mb.plume_mass)\n", + "plt.xlabel(\"time (days)\")\n", + "plt.ylabel(\"plume mass (g)\")\n", + "plt.show()\n", + "plt.plot(mb.t, mb.degraded_mass)\n", + "plt.xlabel(\"time (days)\")\n", + "plt.ylabel(\"degraded mass (g)\")\n", + "plt.title(\"Degraded plume mass of instant reaction model, compared to no decay.\")\n", + "plt.show()" + ] + } + ], + "metadata": { + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 57134ef7353e95b016313aed5d2c5a0999fb96fa Mon Sep 17 00:00:00 2001 From: raar1 Date: Wed, 4 Feb 2026 14:31:11 +0100 Subject: [PATCH 07/10] Resolve gitignore conflict --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index eb877a7..68c7342 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,8 @@ coverage.xml .tox docs/_build +docs/examples/*.gif + site # ide From c67b0606c75a26432f76329cb3bb773e8d799512 Mon Sep 17 00:00:00 2001 From: raar1 Date: Wed, 4 Feb 2026 14:40:15 +0100 Subject: [PATCH 08/10] Add missing comma in pyproject... --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 7ba4311..c40a714 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -84,8 +84,8 @@ doc = [ "mkdocstrings[python] ==0.27.0", "mdx-include >=1.4.1", "mkdocs-markdownextradata-plugin >=0.2.5", - "mike >=2.1.3" - "mkdocs-jupyter >=0.25.1" + "mike >=2.1.3", + "mkdocs-jupyter >=0.25.1", ] [project.urls] From cac4ca9d86ceddc7c2f4a9e763890d3a88fc1f77 Mon Sep 17 00:00:00 2001 From: Jaro Camphuijsen Date: Wed, 4 Feb 2026 15:13:00 +0100 Subject: [PATCH 09/10] Update docs colorscheme to mibirem colors --- docs/assets/custom.css | 7 ++++++- mkdocs.yml | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/docs/assets/custom.css b/docs/assets/custom.css index 2d5f258..faf9bd2 100644 --- a/docs/assets/custom.css +++ b/docs/assets/custom.css @@ -1,4 +1,9 @@ - +[data-md-color-scheme="default"] { + --md-primary-fg-color: #0071bc; +} +[data-md-color-scheme="slate"] { + --md-primary-fg-color: #699812; +} code.highlight { font-size: 18px; } diff --git a/mkdocs.yml b/mkdocs.yml index 9221d03..aa62adf 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -72,7 +72,7 @@ plugins: show_submodules: true - mkdocs-jupyter: execute: true - include_source: True + include_source: true custom_mathjax_url: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/latest.js?config=TeX-AMS_CHTML-full,Safe" # Styled blocks: https://squidfunk.github.io/mkdocs-material/reference/admonitions/#supported-types From cf56df524b0764414cdbbc517cdbf388dd973c6c Mon Sep 17 00:00:00 2001 From: raar1 Date: Wed, 4 Feb 2026 16:16:58 +0100 Subject: [PATCH 10/10] Add a (hacky) CI step to set notebook execution to false during documentation test builds --- .github/workflows/documentation-test.yml | 4 ++++ ...vert_notebook_execute_false_mkdocs_yaml.py | 20 +++++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 __CI_convert_notebook_execute_false_mkdocs_yaml.py diff --git a/.github/workflows/documentation-test.yml b/.github/workflows/documentation-test.yml index ba7df97..ee13bdc 100644 --- a/.github/workflows/documentation-test.yml +++ b/.github/workflows/documentation-test.yml @@ -23,6 +23,10 @@ jobs: - name: Install dependencies run: | pip install ".[doc]" + - name: Set notebook execution to false (for docs testing) + run: | + pip install pyyaml + python __CI_convert_notebook_execute_false_mkdocs_yaml.py mkdocs.yml - name: Test mkdocs build env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/__CI_convert_notebook_execute_false_mkdocs_yaml.py b/__CI_convert_notebook_execute_false_mkdocs_yaml.py new file mode 100644 index 0000000..cfaf1ce --- /dev/null +++ b/__CI_convert_notebook_execute_false_mkdocs_yaml.py @@ -0,0 +1,20 @@ +import sys +import yaml + +if len(sys.argv) != 2: + sys.exit('Usage: python __CI_convert_notebook_execute_false_mkdocs_yaml.py OUTPUTFILE.yml') + +outfname = sys.argv[1] +print(f"Setting mkdocs-jupyter execute flag to false in mkdocs.yml, writing output to {outfname}") + + +with open('mkdocs.yml', 'r') as infile: + conf = yaml.safe_load(infile) + +# If plugins block contains a dict with key 'mkdocs-jupyter', set the 'execute' flag in there to false +nb_dict_list = list(filter(lambda param: 'mkdocs-jupyter' in param, conf['plugins'])) +if len(nb_dict_list) == 1: + nb_dict_list[0]['mkdocs-jupyter']['execute'] = False + + with open(outfname, 'w') as outfile: + yaml.safe_dump(conf, outfile, sort_keys=False)