{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Model Creation\nIn this example, we'll create an iMODFLOW model\nfrom scratch with complex boundary conditions\nand horizontal barriers.\n\nThere are two surface water systems:\nthe outer two edges of the grid feature ditches\nwith a rising stage, whereas the central ditch\nhas a periodic boundary conditions with a\nsummer and winter stage.\n\nThe model will be written as a projectfile\nwith a set of IDFs containing all the grid information\nand a .tim file containing the time discretization.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll start with the usual imports, supplied\nwith geopandas and shapely to specify vector data\nfor the hfb package.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import geopandas as gpd\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\nfrom shapely.geometry import LineString\n\nimport imod\nimport imod.flow as flow" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Discretization\n\nWe'll start off by creating a model discretization.\nThe model consists of a 3 by 9 by 9 three-dimensional grid.\n\nWe'll specify the grid first.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "shape = nlay, nrow, ncol = 3, 9, 9\n\ndx = 100.0\ndy = -100.0\ndz = np.array([5, 30, 100])\nxmin = 0.0\nxmax = dx * ncol\nymin = 0.0\nymax = abs(dy) * nrow\ndims = (\"layer\", \"y\", \"x\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll create the coordinates\nwhich set the grid dimensions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "layer = np.arange(1, nlay + 1)\ny = np.arange(ymax, ymin, dy) + 0.5 * dy\nx = np.arange(xmin, xmax, dx) + 0.5 * dx\ncoords = {\"layer\": layer, \"y\": y, \"x\": x}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The vertical grid discretization (tops and bottoms) are set with a 1D DataArray.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "surface = 0.0\ninterfaces = np.insert((surface - np.cumsum(dz)), 0, surface)\n\nbottom = xr.DataArray(interfaces[1:], coords={\"layer\": layer}, dims=\"layer\")\ntop = xr.DataArray(interfaces[:-1], coords={\"layer\": layer}, dims=\"layer\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll have to create a time discretization as well.\nCreate 1 year of monthly timesteps\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "times = pd.date_range(start=\"1/1/2018\", end=\"12/1/2018\", freq=\"MS\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll create our first 3 dimensional grid here,\nthe ibound grid, which sets where active cells are (ibound = 1.0)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ibound = xr.DataArray(np.ones(shape), coords=coords, dims=dims)\nibound" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hydrogeology\n\nWe'll create a very simple hydrogeology,\nby specifying kh, kva, and sto as constants\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kh = 10.0\nkva = 1.0\nsto = 0.001" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial conditions\n\nWe do not put much effort in the creation of the initial conditions\nin this example, instead we copy the ibounds.\nThis is a 3D grid filled with the value 1, and we can use it as a\ninital condition as well.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "starting_head = ibound.copy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Boundary conditions\n\nWe will put some more effort in creating some complex\nboundary conditions. We'll create both two outer ditches\nwith a rising stage, as well as a central ditch with periodic\n(summer-winter) stage.\n\n### Rising outer ditches\n\nFirst, we'll create rising trend,\nby creating a 1D array ones with the same size as the time dimension,\nand computing the cumulative sum over it.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "trend = np.ones(times[:-1].shape)\ntrend = np.cumsum(trend)\ntrend_da = xr.DataArray(trend, coords={\"time\": times[:-1]}, dims=[\"time\"])\n\ntrend_da" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we assign values only to edges of model x domain.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "is_x_edge = starting_head.x.isin([x[0], x[-1]])\nhead_edge = starting_head.where(is_x_edge)\nhead_edge" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's multiply our\n1D DataArray with dimension time, with the static 3D grid\nwith dimension layer, y, x,\nwhich xarray automatically broadcasts to a 4D array,\nwith dimensions time, layer, y, x\nThis finishes our\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "head_edge_rising = trend_da * head_edge\nhead_edge_rising" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Periodic central ditch\n\nWe'll take only the central column of the grid with\n(where), the rest will be set to np.nan,\nand from this we'll select only the upper layer,\nas the ditch will be located only in the upper layer.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "is_x_central = starting_head.x == x[4]\nhead_central = starting_head.where(is_x_central).sel(layer=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create period times, we let these times start before\nthe model starts.\nThis is necessary because iMODFLOW only forward fills periods.\nOtherwise, in this case there wouldn't be a\nperiodic boundary condition until april.\n\nWe will do this by selecting the months april and october,\nand then subtracting a year\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "period_times = times[[3, 9]] - np.timedelta64(365, \"D\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are creating a summer and winter level.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "periods_da = xr.DataArray([4, 10], coords={\"time\": period_times}, dims=[\"time\"])\nhead_periodic = periods_da * head_central\n\nhead_periodic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create dictionary to tell iMOD\nwhich period name corresponds to which date.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timemap = {\n period_times[0]: \"summer\",\n period_times[1]: \"winter\",\n}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Wells\n\nWells are specified as a pandas dataframe.\nWe create a diagonal line of wells through the domain.\n\nBecause we can.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "wel_df = pd.DataFrame()\nwel_df[\"id_name\"] = np.arange(len(x)).astype(str)\nwel_df[\"x\"] = x\nwel_df[\"y\"] = y\nwel_df[\"rate\"] = dx * dy * -1 * 0.5\nwel_df[\"time\"] = np.tile(times[:-1], 2)[: len(x)]\nwel_df[\"layer\"] = 2\n\nwel_df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Horizontal Flow Barrier\n\nCreate barriers between ditches in layer 1 and 2 (but not 3).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "line1 = LineString([(x[2], ymin), (x[2], ymax)])\nline2 = LineString([(x[7], ymin), (x[7], ymax)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll have to repeat each line for each layer\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "lines = np.array([line1, line2, line1, line2], dtype=\"object\")\nhfb_layers = np.array([1, 1, 2, 2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can specify names for our own bookkeeping\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "id_name = [\"left_upper\", \"right_upper\", \"left_lower\", \"right_lower\"]\n\n# The hfb has to specified as a geopandas GeoDataFrame\n# _\n\nhfb_gdf = gpd.GeoDataFrame(\n geometry=lines, data=dict(id_name=id_name, layer=hfb_layers, resistance=100.0)\n)\n\nhfb_gdf" ] }, { "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 = flow.ImodflowModel(\"my_first_imodflow_model\")\nm[\"pcg\"] = flow.PreconditionedConjugateGradientSolver()\n\nm[\"bnd\"] = flow.Boundary(ibound)\nm[\"top\"] = flow.Top(top)\nm[\"bottom\"] = flow.Bottom(bottom)\n\nm[\"khv\"] = flow.HorizontalHydraulicConductivity(kh)\nm[\"kva\"] = flow.VerticalAnisotropy(kva)\nm[\"sto\"] = flow.StorageCoefficient(sto)\n\nm[\"shd\"] = flow.StartingHead(starting_head)\n\nm[\"chd\"] = flow.ConstantHead(head=head_edge_rising)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create periodic boundary condition\nand specify it as a periodic stress package.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m[\"ghb\"] = flow.GeneralHeadBoundary(head=head_periodic, conductance=10.0)\nm[\"ghb\"].periodic_stress(timemap)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can specify a second stress package as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m[\"ghb2\"] = flow.GeneralHeadBoundary(head=head_periodic + 10.0, conductance=1.0)\n# You also need to specify periodic stresses for second system.\nm[\"ghb2\"].periodic_stress(timemap)\n\nm[\"wel\"] = flow.Well(**wel_df)\n\nm[\"hfb\"] = flow.HorizontalFlowBarrier(**hfb_gdf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Specify output control; -1 specifies to save the output of interest for all layers.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m[\"oc\"] = flow.OutputControl(save_head=-1, save_flux=-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "imod-python wants you to provide the model endtime to your time_discretization!\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m.create_time_discretization(additional_times=times[-1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we write the model\nWrites both .IDFs as well as projectfile, an inifile,\nand a .tim file that contains the time discretization.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "modeldir = imod.util.temporary_directory()\nm.write(directory=modeldir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run\n\nYou can run the model using the comand prompt and the iMOD executables.\nThis is part of the iMOD v5 release, which can be downloaded here:\nhttps://oss.deltares.nl/web/imod/download-imod5 .\niMOD only works on Windows.\n\nTo run your model, open up a command prompt\nand run the following commands:\n\n.. code-block:: batch\n\n cd c:\\path\\to\\modeldir\n c:\\path\\to\\imod\\folder\\iMOD_v5_3_X64R.EXE config_run.ini\n c:\\path\\to\\imod\\folder\\iMODFLOW_V5_3_METASWAP_SVN1977_X64R.EXE my_first_imodflow_model.nam\n\n

#### Note

iMODFLOW requires you to change directory into the model directory\n before calling its executable. Otherwise the program will throw an error.

\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.5" } }, "nbformat": 4, "nbformat_minor": 0 }