{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Regional model\n\nThis example shows a simplified script for building a groundwater model in the\nnortheast of the Netherlands. A primary feature of this area is an ice-pushed\nridge called the Hondsrug. This examples demonstrates modifying external data\nfor use in a MODFLOW6 model.\n\nIn overview, the model features:\n\n    * Thirteen layers: seven aquifers and six aquitards;\n    * A dense ditch network in the east;\n    * Pipe drainage for agriculture;\n    * Precipitation and evapotranspiration.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        ""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll start with the usual imports, and an import from scipy.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport scipy.ndimage\nimport xarray as xr\n\nimport imod"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Before starting to create the input data, we will create the groundwater\nmodel variable (gwf_model) using [imod.mf6.GroundwaterFlowModel](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=groundwater%20flow%20model#imod.mf6.GroundwaterFlowModel).\nThe data from all the model packages will be added to this variable.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model = imod.mf6.GroundwaterFlowModel()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This package allows specifying a regular MODFLOW grid. This grid is assumed\nto be rectangular horizontally, but can be distorted vertically.\n\n## Load data\n\nWe'll load the data from the examples that come with this package.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "layermodel = imod.data.hondsrug_layermodel()\n\n# Make sure that the idomain is provided as integers\nidomain = layermodel[\"idomain\"].astype(int)\n\n# We only need to provide the data for the top as a 2D array. Modflow 6 will\n# compare the top against the uppermost active bottom cell.\ntop = layermodel[\"top\"].max(dim=\"layer\")\n\nbot = layermodel[\"bottom\"]\n\ntop.plot.imshow()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding information to the DIS package\n\nThe following step is to add the previously created discretization data to\nthe gwf_model variable.  This is done using the function\n[imod.mf6.StructuredDiscretization](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.StructuredDiscretization).\nThe data to include is the top of the model domain, the bottom of the layers,\nand the idomain. All this information comes from the previously imported\ntifs (now converted to [xarray.DataArray](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.html#xarray.DataArray).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"dis\"] = imod.mf6.StructuredDiscretization(\n    top=top, bottom=bot, idomain=idomain\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Node property flow package - NPF\n\nThis package contains the information related to the aquifer properties used to calculate\nhydraulic conductance. This package replaces the Layer Property Flow (LPF),\nBlock-Centered Flow (BCF), and Upstream Weighting (UPW) packages from previous MODFLOW versions.\n\n## Hydraulic conductivity\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "k = layermodel[\"k\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## icelltype\n\nThe cell type to be used in the model (confined or convertible) can be\ndefined under ICELLTYPE, which is an input to the NPF package.  ICELLTYPE ==\n0: *Confined cell* - Constant transmissivity ICELLTYPE != 0: *Convertible\ncell* - Transmissivity varies depending on the calculated head in the cell\n(based on the saturated cell thickness)\n\nIn this example, all layers have a ICELLTYPE equal to 0, indicating confined\ncells.  This is defined in the following section.\n\n## Adding information to the NPF package\n\nThe information for the NPF package is added to the gwf_model variable using\n[imod.mf6.NodePropertyFlow](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.NodePropertyFlow).\nThe information included is the icelltype value (equal to zero), the array\nfor  the hydraulic conductivity (considered to be the same for the horizontal\nand vertical direction) and, optionally, the\nvariable_vertical_conductance, dewatered, perched and save_flows options have\nbeen activated.  For more details about the meaning of these variables and\nother variables available to be used within this package, please refer to the\n[documentation](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.NodePropertyFlow).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"npf\"] = imod.mf6.NodePropertyFlow(\n    icelltype=0,\n    k=k,\n    k33=k,\n    variable_vertical_conductance=True,\n    dewatered=True,\n    perched=True,\n    save_flows=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Initial conditions package - IC\n\nThis package reads the starting heads for a simulation.\n\n## Starting heads interpolation\n\nThe starting heads to be used in this model are based on the interpolation of\nx-y head measurements, which were interpolated on a larger area.  This\nexample was created in this example --insert-link-here--\n\nThe heads were interpolated on a larger area, therefore these have to be\nclipped first\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "initial = imod.data.hondsrug_initial()\ninterpolated_head_larger = initial[\"head\"]\n\nxmin = 237_500.0\nxmax = 250_000.0\nymin = 559_000.0\nymax = 564_000.0\n\ninterpolated_head = interpolated_head_larger.sel(\n    x=slice(xmin, xmax), y=slice(ymax, ymin)\n)\n\n# Plotting the clipped interpolation\nfig, ax = plt.subplots()\ninterpolated_head.plot.imshow(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The final step is to assign the 2D heads interpolation to all the\nmodel layers (as a reference value) by using the xarray tool\n[xarray.full_like](http://xarray.pydata.org/en/stable/generated/xarray.full_like.html#xarray.full_like).\nThe 3d idomain array is used as reference for the geometry and then\nits original values are replaced by NaNs.\nThis array is combined with the interpolated_head array using the xarray\n[DataArray.combine_first](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.combine_first.html#xarray.DataArray.combine_first)\noption.\nThe final result is an starting_heads xarray where all layers have the 2d interpolated_head information.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Assign interpolated head values to all the model layers\nlike_3d = xr.full_like(idomain, np.nan, dtype=float)\nstarting_head = like_3d.combine_first(interpolated_head)\n# Consequently ensure no data is specified in inactive cells:\nstarting_head = starting_head.where(idomain == 1)\n\nstarting_head"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding information to the IC package\n\nThe function for indicating the initial conditions is\n[imod.mf6.InitialConditions](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.InitialConditions).\nIt is necessary to indicate the value(s) to be considered as the initial\n(starting) head of the simulation.\nIn this case, this value is equal to the previously created starting_head array.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"ic\"] = imod.mf6.InitialConditions(starting_head)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Constant head package - CHD\n\nThis package allows to indicate if the head varies with time,\nif it is constant or if it is inactive.\n\n## Constant head edge\n\nThe previously interpolated starting_head array will be used to define\nthe constant head value which will be used along the model boundaries.\nA function is defined to indicate the location of the outer edge\n(returning a boolean array).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def outer_edge(da):\n    data = da.copy()\n    from_edge = scipy.ndimage.binary_erosion(data)\n    is_edge = (data == 1) & (from_edge == 0)\n    return is_edge.astype(bool)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the next calculations, it is necessary to create a template array\nwhich can be used for assigning the corresponding geometry to other arrays.\nIn this case, a 2d template is created based on the idomain layer information\nand filled with ones.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "like_2d = xr.full_like(idomain.isel(layer=0), 1)\nlike_2d"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using the previously created function and the 2d template,\nthe outer edge is defined for this example.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "edge = outer_edge(xr.full_like(like_2d.drop_vars(\"layer\"), 1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding information to the CHD package\n\nTo add the information to the CHD package within the gwf_model variable, the\n[imod.mf6.ConstantHead](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.ConstantHead)\nfunction is used.\nThe required information is the head array for this boundary condition.\nIn this example, the starting_head array is selected where the idomain is > 0 (active)\nand it is located in the edge array.\n\nIt is also possible (and optional) to indicate if the CHD information will be written\nto the listing file after it is read (print_input), if the constant head flow rates will\nbe printed to the listing file for every stress period time step\nin which \u201cBUDGET PRINT\u201d is specified in Output Control (print_flows)\nand if the constant head flow terms will be written to the file\nspecified with \u201cBUDGET FILEOUT\u201d in Output Control (save_flows).\nBy default, these three options are set to False.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"chd\"] = imod.mf6.ConstantHead(\n    starting_head.where((idomain > 0) & edge),\n    print_input=False,\n    print_flows=True,\n    save_flows=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Recharge\n\nThis package is used to represent areally distributed recharge to the groundwater system.\nTo calculate the recharge, the precipitation and evapotranspiration\ninformation from the KNMI website has been downloaded for the study area.\nThis information is saved in netCDF files, which have been imported using the\nxarray function\n[xr.open_dataset](http://xarray.pydata.org/en/stable/generated/xarray.open_dataset.html#xarray.open_dataset),\nslicing the area to the model's miminum and maximum dimensions.\n\nNote that the meteorological data has mm/d as unit and\nthis has to be converted to m/d for Modflow 6.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xmin = 230_000.0\nxmax = 257_000.0\nymin = 550_000.0\nymax = 567_000.0\n\nmeteorology = imod.data.hondsrug_meteorology()\npp = meteorology[\"precipitation\"]\nevt = meteorology[\"evapotranspiration\"]\n\npp = pp.sel(x=slice(xmin, xmax), y=slice(ymax, ymin)) / 1000.0  # from mm/d to m/d\nevt = evt.sel(x=slice(xmin, xmax), y=slice(ymax, ymin)) / 1000.0  # from mm/d to m/d"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Recharge - Steady state\n\nFor the steady state conditions of the model,\nthe data from the period 2000 to 2009 was considered as reference.\nThe initial information was sliced to this time period and averaged\nto obtain the a mean value grid. This process was done for both\nprecipitation and evapotranspiration datasets.\n\n**Precipitation**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pp_ss = pp.sel(time=slice(\"2000-01-01\", \"2009-12-31\"))\npp_ss_mean = pp_ss.mean(dim=\"time\")\n\nfig, ax = plt.subplots()\npp_ss_mean.plot(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Evapotranspiration**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "evt_ss = evt.sel(time=slice(\"2000-01-01\", \"2009-12-31\"))\nevt_ss_mean = evt_ss.mean(dim=\"time\")\n\nfig, ax = plt.subplots()\nevt_ss_mean.plot(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the recharge calculation, a first estimate\nis the difference between the precipitation and evapotranspiration values.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_ss = pp_ss_mean - evt_ss_mean\n\nfig, ax = plt.subplots()\nrch_ss.plot.imshow(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Recharge - Transient\n\nThe transient model will encompass the period from 2010 to 2015.\nThe initial pp and evt datasets have been sliced to this time frame.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pp_trans = pp.sel(time=slice(\"2010-01-01\", \"2015-12-31\"))\nevt_trans = evt.sel(time=slice(\"2010-01-01\", \"2015-12-31\"))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As previously done, it is assumed that the recharge is equal\nto the difference between precipitation and evapotranspiration as a first estimate.\nFurthermore, the negative values found after doing this calculation have been\nreplaced by zeros, as the recharge should not have a negative value.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_trans = pp_trans - evt_trans\nrch_trans = rch_trans.where(rch_trans > 0, 0)  # check negative values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The original information is on a daily step, so it is going to be\nresampled to a yearly step by using the xarray function\n[Dataset.resample](http://xarray.pydata.org/en/stable/generated/xarray.Dataset.resample.html#xarray.Dataset.resample).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_trans_yr = rch_trans.resample(time=\"A\", label=\"left\").mean()\nrch_trans_yr"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To create the final recharge for the transient simulation,\nthe steady state information needs to be concatenated to the transient recharge data.\nThe steady state simulation will be run for one second.\nThis is achieved by using the numpy\n[Timedelta function](https://numpy.org/doc/stable/reference/arrays.datetime.html),\nfirst creating a time delta of 1 second, which is assigned to the steady state recharge information.\nThis dataset is then concatenated using the xarray function\n[xarray.concat](http://xarray.pydata.org/en/stable/generated/xarray.concat.html#xarray.concat)\nto the transient information and indicating that the dimension to join is \"time\".\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "starttime = \"2009-12-31\"\n\n# Add first steady-state\ntimedelta = np.timedelta64(1, \"s\")  # 1 second duration for initial steady-state\nstarttime_steady = np.datetime64(starttime) - timedelta\nrch_ss = rch_ss.assign_coords(time=starttime_steady)\n\nrch_ss_trans = xr.concat([rch_ss, rch_trans_yr], dim=\"time\")\nrch_ss_trans"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The data obtained from KNMI has different grid dimensions\nthan the one considered in this example. To fix this,\nimod-python includes the option\n[imod.prepare.Regridder](https://deltares.gitlab.io/imod/imod-python/api/prepare.html#imod.prepare.Regridder),\nwhich modifies the original grid dimensions to a different one.\nIt is also possible to define the regridding method such as\n``nearest``, ``multilinear``, ``mean``, among others.\nIn this case, ``mean`` was selected and the 2d template (like_2d)\nwas used as reference, as this is the geometry to be considered in the model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_ss_trans = imod.prepare.Regridder(method=\"mean\").regrid(rch_ss_trans, like_2d)\nrch_ss_trans"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The previously created recharge array is a 2D array\nthat needs to be assigned to a 3D array. This is done using the xarray\n[DataArray.where](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.where.html#xarray.DataArray.where)\noption, where the recharge values are applied to the cells where the\nidomain value is larger than zero (that is, the active cells) and for the uppermost\nactive cell (indicated by the minimum layer number).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_total = rch_ss_trans.where(\n    idomain[\"layer\"] == idomain[\"layer\"].where(idomain > 0).min(\"layer\")\n)\nrch_total"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, transposing the array dimensions using\n[DataArray.transpose](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.transpose.html#xarray.DataArray.transpose)\nso they are in the correct order.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_total = rch_total.transpose(\"time\", \"layer\", \"y\", \"x\")\nrch_total\n\nfig, ax = plt.subplots()\nrch_total.isel(layer=2, time=6).plot.imshow(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding information to the RCH package\n\nThe information for the RCH package is added with the function\n[imod.mf6.Recharge](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.Recharge).\nIt is required to insert the recharge flux rate, and it is optional\nto include the print_input, print_flows and save_flows information.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"rch\"] = imod.mf6.Recharge(rch_total)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Drainage package - DRN\n\nThe drain package is used to simulate features that remove water from the aquifer,\nsuch as agricultural drains or springs.\nThis occurs at a rate proportional to the head difference between the head in the\naquifer and the drain elevation\n(the head in the aquifer has to be above that elevation).\nThe conductance is the proportionality constant.\n\n## Import drainage information\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "drainage = imod.data.hondsrug_drainage()\n\npipe_cond = drainage[\"conductance\"]\npipe_elev = drainage[\"elevation\"]\n\npipe_cond"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding information to the DRN package\n\nTo add the information to the DRN package within the gwf_model variable, the\n[imod.mf6.Drainage](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.Drainage)\nfunction is used. It is required to add the previously created arrays for\nthe drain elevation and the drain conductance.\nIt is optional to insert the information for\n``print_input``, ``print_flows`` and ``save_flows``\nwhich are set to False by default.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"drn-pipe\"] = imod.mf6.Drainage(conductance=pipe_cond, elevation=pipe_elev)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# River package - RIV\n\nThis package simulates the effects of flow between\nsurface-water features and groundwater systems.\n\n## Import river information\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "river = imod.data.hondsrug_river()\nriv_cond = river[\"conductance\"]\nriv_stage = river[\"stage\"]\nriv_bot = river[\"bottom\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding information to the RIV package\n\nThe data is assigned to the gwf_model variable by using\n[imod.mf6.River](https://deltares.gitlab.io/imod/imod-python/api/mf6.html#imod.mf6.River),\nbased on the previously imported conductance, stage and bottom arrays.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"riv\"] = imod.mf6.River(\n    conductance=riv_cond, stage=riv_stage, bottom_elevation=riv_bot\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Storage package - STO\n\nWhen the STO Package is included in a model, storage changes\nwill be calculated, and thus, the model will be transient.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ss = 0.0003\nlayer = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])\nsy = xr.DataArray(\n    [0.16, 0.16, 0.16, 0.16, 0.15, 0.15, 0.15, 0.15, 0.14, 0.14, 0.14, 0.14, 0.14],\n    {\"layer\": layer},\n    (\"layer\",),\n)\ntimes_sto = np.array(\n    [\n        \"2009-12-30T23:59:59.00\",\n        \"2009-12-31T00:00:00.00\",\n        \"2010-12-31T00:00:00.00\",\n        \"2011-12-31T00:00:00.00\",\n        \"2012-12-31T00:00:00.00\",\n        \"2013-12-31T00:00:00.00\",\n        \"2014-12-31T00:00:00.00\",\n    ],\n    dtype=\"datetime64[ns]\",\n)\n\ntransient = xr.DataArray(\n    [False, True, True, True, True, True, True], {\"time\": times_sto}, (\"time\",)\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding information to the STO package\n\nThe data is assigned to the gwf_model variable by using\n[imod.mf6.Storage](https://deltares.gitlab.io/imod/imod-python/api/generated/mf6/imod.mf6.SpecificStorage.html).\nIt is necessary to indicate the values of specific storage,\nspecific yield and if the layers are convertible.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"sto\"] = imod.mf6.SpecificStorage(\n    specific_storage=ss,\n    specific_yield=sy,\n    transient=transient,\n    convertible=0,\n    save_flows=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Output Control package - OC\n\nThis package determines how and when heads are printed to the\nlisting file and/or written to a separate binary output file\n\n## Adding information to the OC package\n\nThe function\n[imod.mf6.OutputControl](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.OutputControl)\nis used to store the information for this package.\nIt is possible to indicate if the heads and budget information is saved\nat the end of each stress period (``last``),\nfor all timesteps a stress period (``all``),\nor at the start of a stress period (``first``)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model[\"oc\"] = imod.mf6.OutputControl(save_head=\"last\", save_budget=\"last\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Model simulation\n\nIn MODFLOW 6, the concept of \"model\" is that part of the program\nthat solves a hydrologic process. MODFLOW 6 documentation supports\none type of model: the GWF Model.\nIt is possible within the MODFLOW 6 framewotk to solve multiple,\ntightly coupled, numerical models in a single system of equation,\nwhich may be multiple models of the same type or of different types.\n\nThe previously created gwf_model variable now contains\nthe information from all the variables.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Attach the model information to a simulation\n\nThe function\n[imod.mf6.Modflow6Simulation](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.Modflow6Simulation)\nallows to assign models to a simulation (in this case, the gwf_model).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simulation = imod.mf6.Modflow6Simulation(\"mf6-mipwa2-example\")\nsimulation[\"GWF_1\"] = gwf_model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solver settings\n\nThe solver settings are indicated using\n[imod.mf6.Solution](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=structured%20discretization#imod.mf6.Solution).\nIf the values are not indicated manually, the defaults values will be considered.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simulation[\"solver\"] = imod.mf6.Solution(\n    modelnames=[\"GWF_1\"],\n    print_option=\"summary\",\n    csv_output=False,\n    no_ptc=True,\n    outer_dvclose=1.0e-4,\n    outer_maximum=500,\n    under_relaxation=None,\n    inner_dvclose=1.0e-4,\n    inner_rclose=0.001,\n    inner_maximum=100,\n    linear_acceleration=\"cg\",\n    scaling_method=None,\n    reordering_method=None,\n    relaxation_factor=0.97,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Assign time discretization\n\nThe time discretization of this model is 6 years.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simulation.create_time_discretization(\n    additional_times=[\"2009-12-30T23:59:59.000000000\", \"2015-12-31T00:00:00.000000000\"]\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run the model\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The following lines assume the ``mf6`` executable is available on your PATH.\n  `The Modflow 6 examples introduction <mf6-introduction>` shortly\n  describes how to add it to yours.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "modeldir = imod.util.temporary_directory()\nsimulation.write(modeldir, binary=False)\nsimulation.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Results visualization\n\nThe next section indicated how to visualize the model results.\n\n## Import heads results\n\nThe heads results are imported using\n[imod.mf6.open_hds](https://deltares.gitlab.io/imod/imod-python/api/mf6.html?highlight=imod%20mf6%20open_hds#imod.mf6.open_hds)\non the background.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "hds = simulation.open_head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can plot the data of an individual layer as follows\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nhds.sel(layer=3).isel(time=3).plot(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As you can see layer 3 has some missing cells in the west\nWhereas layer 4 only contains active cells in the\neastern peatland area\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nhds.sel(layer=4).isel(time=3).plot(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Layer 5 contains more data towards the west,\nbut has no active cells in the centre.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nhds.sel(layer=5).isel(time=3).plot(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As you can see the data is individual layers\nhave lots of inactive in different places.\n\nIt is difficult for this model to get a good idea\nwhat is happening across the area based on 1 layer alone.\nLuckily xarray allows us to compute the mean across a selection\nof layers and plot this.\n\nBy first selecting 3 layers with ``sel```,\nand then computing the mean across the layer dimension\nwith ``mean(dim=\"layer\")``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nhds.sel(layer=slice(3, 5)).mean(dim=\"layer\").isel(time=3).plot(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Assign dates to head\n\nMODFLOW6 has no concept of a calendar, so the output is not labelled only\nin terms of \"time since start\" in floating point numbers. For this model\nthe time unit is days and we can assign a date coordinate as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "starttime = pd.to_datetime(\"2000-01-01\")\ntimedelta = pd.to_timedelta(hds[\"time\"], \"D\")\nhds = hds.assign_coords(time=starttime + timedelta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Extract head at points\n\nA typical operation is to extract simulated heads at point locations to\ncompare them with measurements. In this example, we select the heads at\ntwo points:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x = [240_000.0, 244_000.0]\ny = [560_000.0, 562_000.0]\nselection = imod.select.points_values(hds, x=x, y=y)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The result can be converted into a pandas dataframe for timeseries analysis,\nor written to a variety of tabular file formats.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dataframe = selection.to_dataframe().reset_index()\ndataframe = dataframe.rename(columns={\"index\": \"id\"})\ndataframe"
      ]
    }
  ],
  "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.11.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}