{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Henry\n\nThis example illustrates how to setup a variable density groundwater flow and\ntransport model using the ``imod`` package and associated packages.\n\nIn overview, we'll set the following steps:\n\n    * Create a suitable 2d (x, z) grid.\n    * Create a groundwater flow model, with variable density.\n    * Create a solute transport model.\n    * Combine these models into a single MODFLOW6 simulation.\n    * Write to modflow6 files.\n    * Run the model.\n    * Open the results back into xarray DataArrays.\n    * Visualize the results.\n\nWe are simulating the Henry problem, although not the original one but the one\noutlined in the MODFLOW 6 manual (jupyter notebooks example 51). This is the\nmodified Henry problem with half the inflow rate.\n\nThe domain is a vertically oriented two dimensional rectangle, which is 2 m\nlong and 1 m high. Water flows in over the left boundary with a fixed rate,\nwhich is represented by a Well package. The right boundary is in direct contact\nwith hydrostatic seawater with a density of 1025 kg m:sup:`-3`. This is\nrepresented by a General Head Boundary package.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        ""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll start with the usual imports. As this is a simple (synthetic)\nstructured model, we can make due with few packages.\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": [
        "We'll start by defining the (vertical) rectangular domain and the physical\nparameters of the model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nlay = 40\nnrow = 1\nncol = 80\nshape = (nlay, nrow, ncol)\n\ntotal_flux = 5.7024  # m3/d\nk = 864.0  # m/d\nporosity = 0.35\nmax_concentration = 35.0\nmin_concentration = 0.0\nmax_density = 1025.0\nmin_density = 1000.0\ndiffusion_coefficient = 0.57024\nlongitudinal_horizontal = 0.1\ntransversal_horizontal1 = 0.01\n\n# Time\nstart_date = pd.to_datetime(\"2020-01-01\")\nduration = pd.to_timedelta(\"0.5d\")\n\n# Domain size\nxmax = 2.0\nxmin = 0.0\ndx = (xmax - xmin) / ncol\nzmin = 0.0\nzmax = 1.0\ndz = (zmax - zmin) / nlay\n\nx = np.arange(xmin, xmax, dx) + 0.5 * dx\ny = np.array([0.5])\nlayer = np.arange(1, 41, 1)\n\ndy = -1.0\ncoords = {\"layer\": layer, \"y\": y, \"x\": x, \"dy\": dy, \"dx\": dx}\ndims = (\"layer\", \"y\", \"x\")\nidomain = xr.DataArray(np.ones(shape, dtype=int), coords=coords, dims=dims)\n\ntop = 1.0\nbottom = xr.DataArray(\n    np.arange(zmin, zmax, dz)[::-1],\n    {\"layer\": layer},\n    (\"layer\",),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now make the flow model. We'll start with the non-boundary condition\npackages.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwf_model = imod.mf6.GroundwaterFlowModel()\ngwf_model[\"dis\"] = imod.mf6.StructuredDiscretization(\n    top=top, bottom=bottom, idomain=idomain\n)\ngwf_model[\"npf\"] = imod.mf6.NodePropertyFlow(\n    icelltype=0,\n    k=k,\n)\ngwf_model[\"sto\"] = imod.mf6.SpecificStorage(\n    specific_storage=1.0e-4,\n    specific_yield=0.15,\n    transient=False,\n    convertible=0,\n)\ngwf_model[\"ic\"] = imod.mf6.InitialConditions(start=0.0)\ngwf_model[\"oc\"] = imod.mf6.OutputControl(save_head=\"last\", save_budget=\"last\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's make the constant head boundary condition.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ghb_head = xr.ones_like(idomain, dtype=float)\nghb_head[:, :, :-1] = np.nan\n\nghb_conc = xr.full_like(idomain, max_concentration, dtype=float)\nghb_conc[:, :, :-1] = np.nan\nghb_conc = ghb_conc.expand_dims(species=[\"salinity\"])\n\nconductance = xr.full_like(idomain, 864.0 * 2.0, dtype=float)\nconductance[:, :, :-1] = np.nan\n\ngwf_model[\"right_boundary\"] = imod.mf6.GeneralHeadBoundary(\n    head=ghb_head,\n    conductance=conductance,\n    concentration=ghb_conc,\n    concentration_boundary_type=\"AUX\",\n    print_input=True,\n    print_flows=True,\n    save_flows=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "... and the constant flux condition.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "flux_concentration = xr.DataArray(\n    data=np.full((1, nlay), min_concentration),\n    dims=[\"species\", \"cell\"],\n    coords={\"species\": [\"salinity\"], \"cell\": layer},\n)\n\ngwf_model[\"left_boundary\"] = imod.mf6.WellDisStructured(\n    layer=layer,\n    row=np.full_like(layer, 1, dtype=int),\n    column=np.full_like(layer, 1, dtype=int),\n    rate=np.full_like(layer, 0.5 * (total_flux / nlay), dtype=float),\n    concentration=flux_concentration,\n    concentration_boundary_type=\"AUX\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since we are solving a variable density problem, we need to add the buoyancy\npackage. It will use the species \"salinity\" that we are simulating with a\ntransport model defined below.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "slope = (max_density - min_density) / (max_concentration - min_concentration)\ngwf_model[\"buoyancy\"] = imod.mf6.Buoyancy(\n    reference_density=min_density,\n    modelname=[\"transport\"],\n    reference_concentration=[min_concentration],\n    density_concentration_slope=[slope],\n    species=[\"salinity\"],\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's make the transport model. It contains the standard packages of\nstorage, dispersion and advection as well as initial condiations and output\ncontrol. Sinks and sources are automatically determined based on packages\nprovided in the flow model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gwt_model = imod.mf6.GroundwaterTransportModel()\ngwt_model[\"ssm\"] = imod.mf6.SourceSinkMixing.from_flow_model(gwf_model, \"salinity\")\ngwt_model[\"adv\"] = imod.mf6.AdvectionTVD()\ngwt_model[\"dsp\"] = imod.mf6.Dispersion(\n    diffusion_coefficient=0.57024,\n    longitudinal_horizontal=0.1,\n    transversal_horizontal1=0.01,\n    xt3d_off=False,\n    xt3d_rhs=False,\n)\n\ngwt_model[\"mst\"] = imod.mf6.MobileStorageTransfer(\n    porosity=porosity,\n)\ngwt_model[\"ic\"] = imod.mf6.InitialConditions(start=max_concentration)\ngwt_model[\"oc\"] = imod.mf6.OutputControl(save_concentration=\"last\", save_budget=\"last\")\ngwt_model[\"dis\"] = gwf_model[\"dis\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "now let's define a simulation using the flow and transport models.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Attach it to a simulation\nsimulation = imod.mf6.Modflow6Simulation(\"henry\")\n\nsimulation[\"flow\"] = gwf_model\nsimulation[\"transport\"] = gwt_model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define solver settings. We need to define separate solutions for the flow and\ntransport models. In this case, we'll use the same settings, but generally\nconvergence settings should differ: the transport model has very different\nunits from the flow model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simulation[\"flow_solver\"] = imod.mf6.Solution(\n    modelnames=[\"flow\"],\n    print_option=\"summary\",\n    csv_output=False,\n    no_ptc=True,\n    outer_dvclose=1.0e-6,\n    outer_maximum=500,\n    under_relaxation=None,\n    inner_dvclose=1.0e-6,\n    inner_rclose=1.0e-10,\n    inner_maximum=100,\n    linear_acceleration=\"bicgstab\",\n    scaling_method=None,\n    reordering_method=None,\n    relaxation_factor=0.97,\n)\nsimulation[\"transport_solver\"] = imod.mf6.Solution(\n    modelnames=[\"transport\"],\n    print_option=\"summary\",\n    csv_output=False,\n    no_ptc=True,\n    outer_dvclose=1.0e-6,\n    outer_maximum=500,\n    under_relaxation=None,\n    inner_dvclose=1.0e-6,\n    inner_rclose=1.0e-10,\n    inner_maximum=100,\n    linear_acceleration=\"bicgstab\",\n    scaling_method=None,\n    reordering_method=None,\n    relaxation_factor=0.97,\n)\n# Collect time discretization\ntimes = [start_date, start_date + duration]\nsimulation.create_time_discretization(additional_times=times)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Increase the number of time steps for the single stress period:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simulation[\"time_discretization\"][\"n_timesteps\"] = 500"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll create a new directory in which we will write and run the model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "modeldir = imod.util.temporary_directory()\nsimulation.write(modeldir, binary=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run the model\n\nThis takes about 20 seconds.\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": [
        "simulation.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Open the results\n\nWe'll open the head and concentration files.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "head = imod.mf6.out.open_hds(\n    modeldir / \"flow/flow.hds\",\n    modeldir / \"flow/dis.dis.grb\",\n)\nconc = imod.mf6.out.open_conc(\n    modeldir / \"transport/transport.ucn\",\n    modeldir / \"flow/dis.dis.grb\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualize the results\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "head.isel(y=0, time=-1).plot.contourf(yincrease=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can check the concentration to see that a fresh-saline interface has been\nformed:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "conc.isel(y=0, time=-1).plot.contourf(yincrease=False, cmap=\"RdYlBu_r\")"
      ]
    }
  ],
  "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.11"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}