{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Hydrocoin\n\nA 2D case from the Hydrological Code Intercomparison (Hydrocoin).\n\nFor more information see:\n\nKonikow, L. F., Sanford, W. E., & Campbell, P. J. (1997).\nConstant-concentration boundary condition: Lessons from the HYDROCOIN\nvariable-density groundwater benchmark problem. Water Resources Research, 33\n(10), 2253-2261. https://doi.org/10.1029/97WR01926\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll start with the usual imports\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport pandas as pd\nimport xarray as xr\n\nimport imod" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Discretization\n\nWe'll start off by creating a model discretization, since\nthis is a simple conceptual model.\nThe model is a 2D cross-section, hence nrow = 1.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nrow = 1 # number of rows\nncol = 45 # number of columns\nnlay = 76 # number of layers\n\ndz = 4.0\ndx = 20.0\ndy = -dx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up tops and bottoms\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "top1D = xr.DataArray(\n np.arange(nlay * dz, 0.0, -dz), {\"layer\": np.arange(1, nlay + 1)}, (\"layer\")\n)\n\nbot = top1D - dz" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up ibound, which sets where active cells are (ibound = 1.0).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bnd = xr.DataArray(\n data=np.full((nlay, nrow, ncol), 1.0),\n coords={\n \"y\": [0.5],\n \"x\": np.arange(0.5 * dx, dx * ncol, dx),\n \"layer\": np.arange(1, 1 + nlay),\n \"dx\": dx,\n \"dy\": dy,\n },\n dims=(\"layer\", \"y\", \"x\"),\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set inactive cells by specifying bnd[index] = 0.0\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bnd[75, :, 0:15] = 0.0\nbnd[75, :, 30:45] = 0.0\n\nfig, ax = plt.subplots()\nbnd.plot(y=\"layer\", yincrease=False, ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Boundary Conditions\n\nSet the constant heads by specifying a negative value in iboud,\nthat is: bnd[index] = -1\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bnd[0, :, :] = -1\n\n\nfig, ax = plt.subplots()\nbnd.plot(y=\"layer\", yincrease=False, ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define WEL data, need to define the x, y, and pumping rate (q)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "weldata = pd.DataFrame()\nweldata[\"x\"] = np.full(1, 0.5 * dx)\nweldata[\"y\"] = np.full(1, 0.5)\nweldata[\"q\"] = 0.28512 # positive, so it's an injection well" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the icbund, which sets which cells\nin the solute transport model are active, inactive or constant.\n\nIn this case the central 15 cells on the top row have a constant concentration,\nAnd, on both sides, the outer 15 cells of the top row are inactive in the transport model.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "icbund = xr.DataArray(\n data=np.full((nlay, nrow, ncol), 1.0),\n coords={\n \"y\": [0.5],\n \"x\": np.arange(0.5 * dx, dx * ncol, dx),\n \"layer\": np.arange(1, nlay + 1),\n \"dx\": dx,\n \"dy\": dy,\n },\n dims=(\"layer\", \"y\", \"x\"),\n)\n\nicbund[75, :, 0:15] = 0.0\nicbund[75, :, 30:45] = 0.0\nicbund[75, :, 15:30] = -1.0\n\nfig, ax = plt.subplots()\nicbund.plot(y=\"layer\", yincrease=False, ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial conditions\n\nDefine the starting concentrations\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sconc = xr.DataArray(\n data=np.full((nlay, nrow, ncol), 0.0),\n coords={\n \"y\": [0.5],\n \"x\": np.arange(0.5 * dx, dx * ncol, dx),\n \"layer\": np.arange(1, nlay + 1),\n \"dx\": dx,\n \"dy\": dy,\n },\n dims=(\"layer\", \"y\", \"x\"),\n)\n\nsconc[75, :, 15:30] = 280.0\n\nfig, ax = plt.subplots()\nsconc.plot(y=\"layer\", yincrease=False, ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define starting heads, these will be inserted in the Basic Flow (BAS) package\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "shd = xr.DataArray(\n data=np.full((nlay, nrow, ncol), 0.0),\n coords={\n \"y\": [0.5],\n \"x\": np.arange(0.5 * dx, dx * ncol, dx),\n \"layer\": np.arange(1, nlay + 1),\n \"dx\": dx,\n \"dy\": dy,\n },\n dims=(\"layer\", \"y\", \"x\"),\n)\n\nshd[0, :, :] = np.array(\n [\n 10,\n 9.772727273,\n 9.545454545,\n 9.318181818,\n 9.090909091,\n 8.863636364,\n 8.636363636,\n 8.409090909,\n 8.181818182,\n 7.954545455,\n 7.727272727,\n 7.5,\n 7.272727273,\n 7.045454545,\n 6.818181818,\n 6.590909091,\n 6.363636364,\n 6.136363636,\n 5.909090909,\n 5.681818182,\n 5.454545455,\n 5.227272727,\n 5,\n 4.772727273,\n 4.545454545,\n 4.318181818,\n 4.090909091,\n 3.863636364,\n 3.636363636,\n 3.409090909,\n 3.181818182,\n 2.954545455,\n 2.727272727,\n 2.5,\n 2.272727273,\n 2.045454545,\n 1.818181818,\n 1.590909091,\n 1.363636364,\n 1.136363636,\n 0.909090909,\n 0.681818182,\n 0.454545455,\n 0.227272727,\n 0.00,\n ]\n)\n\nfig, ax = plt.subplots()\nshd.plot(y=\"layer\", yincrease=False, ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hydrogeology\n\nDefine horizontal hydraulic conductivity\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "khv = xr.DataArray(\n data=np.full((nlay, nrow, ncol), 0.847584),\n coords={\n \"y\": [0.5],\n \"x\": np.arange(0.5 * dx, dx * ncol, dx),\n \"layer\": np.arange(1, nlay + 1),\n \"dx\": dx,\n \"dy\": dy,\n },\n dims=(\"layer\", \"y\", \"x\"),\n)\n\nkhv[75, :, 15:30] = 0.0008475\n\nfig, ax = plt.subplots()\nkhv.plot(y=\"layer\", yincrease=False, ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build\n\nFinally, we build the model.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m = imod.wq.SeawatModel(\"Hydrocoin\")\nm[\"bas\"] = imod.wq.BasicFlow(ibound=bnd, top=304.0, bottom=bot, starting_head=shd)\nm[\"lpf\"] = imod.wq.LayerPropertyFlow(\n k_horizontal=khv, k_vertical=khv, specific_storage=0.0\n)\nm[\"btn\"] = imod.wq.BasicTransport(\n icbund=icbund, starting_concentration=sconc, porosity=0.2\n)\nm[\"adv\"] = imod.wq.AdvectionTVD(courant=1.0)\nm[\"dsp\"] = imod.wq.Dispersion(longitudinal=20.0, diffusion_coefficient=0.0)\nm[\"vdf\"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)\nm[\"wel\"] = imod.wq.Well(\n id_name=\"wel\", x=weldata[\"x\"], y=weldata[\"y\"], rate=weldata[\"q\"]\n)\nm[\"pcg\"] = imod.wq.PreconditionedConjugateGradientSolver(\n max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0\n)\nm[\"gcg\"] = imod.wq.GeneralizedConjugateGradientSolver(\n max_iter=150,\n inner_iter=30,\n cclose=1.0e-6,\n preconditioner=\"mic\",\n lump_dispersion=True,\n)\nm[\"oc\"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)\nm.create_time_discretization(additional_times=[\"2000-01-01T00:00\", \"2010-01-01T00:00\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we write the model, including runfile:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "modeldir = imod.util.temporary_directory()\nm.write(modeldir, resultdir_is_workdir=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run\n\nYou can run the model using the comand prompt and the iMOD-WQ executable.\nThis is part of the iMOD v5 release, which can be downloaded here:\nhttps://oss.deltares.nl/web/imod/download-imod5 .\nThis only works on Windows.\n\nTo run your model, open up a command prompt\nand run the following commands:\n\nbatch\ncd c:\\path\\to\\modeldir\nc:\\path\\to\\imod\\folder\\iMOD-WQ_V5_3_SVN359_X64R.exe Hydrocoin.run\n\nNote that the version name of your executable might differ.\n\n%%\nVisualise results\n-----------------\n\nAfter succesfully running the model, you can\nplot results as follows:\n\n.. code:: python\n\n head = imod.idf.open(modeldir / \"results/head/*.idf\")\n\n fig, ax = plt.subplots()\n head.plot(yincrease=False, ax=ax)\n\n conc = imod.idf.open(modeldir / \"results/conc/*.idf\")\n\n fig, ax = plt.subplots()\n conc.plot(levels=range(0, 35, 5), yincrease=False, ax=ax)\n\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.10.8" } }, "nbformat": 4, "nbformat_minor": 0 }