.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples\imod-wq\Henry-wq.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_imod-wq_Henry-wq.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_imod-wq_Henry-wq.py:


Henry
=====

The classic 2D Henry problem demonstrates the development of a fresh-salt
interface.

.. GENERATED FROM PYTHON SOURCE LINES 10-11

We'll start with the usual imports

.. GENERATED FROM PYTHON SOURCE LINES 11-17

.. code-block:: Python

    import numpy as np
    import pandas as pd
    import xarray as xr

    import imod








.. GENERATED FROM PYTHON SOURCE LINES 18-24

Discretization
--------------

We'll start off by creating a model discretization, since
this is a simple conceptual model.
The model is a 2D cross-section, hence ``nrow = 1``.

.. GENERATED FROM PYTHON SOURCE LINES 24-38

.. code-block:: Python


    nrow = 1
    ncol = 100
    nlay = 50

    dz = 1.0
    dx = 1.0
    dy = -dx

    top1D = xr.DataArray(
        np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, 1 + nlay)}, ("layer")
    )
    bot = top1D - 1.0








.. GENERATED FROM PYTHON SOURCE LINES 39-40

Set up ibound, which sets where active cells are `(ibound = 1.0)`

.. GENERATED FROM PYTHON SOURCE LINES 40-58

.. code-block:: Python

    bnd = xr.DataArray(
        data=np.full((nlay, nrow, ncol), 1.0),
        coords={
            "y": [0.5],
            "x": np.arange(0.5 * dx, dx * ncol, dx),
            "layer": np.arange(1, 1 + nlay),
            "dx": dx,
            "dy": dy,
        },
        dims=("layer", "y", "x"),
    )

    # Boundary Conditions
    # -------------------
    #
    # We define constant head here, after generating the tops, or we'd end up with negative top values
    bnd[:, :, -1] = -1








.. GENERATED FROM PYTHON SOURCE LINES 59-62

Create WEL package

First we scale the pumping rate with discretization

.. GENERATED FROM PYTHON SOURCE LINES 62-64

.. code-block:: Python

    qscaled = 0.03 * (dz * abs(dy))








.. GENERATED FROM PYTHON SOURCE LINES 65-68

Fresh water injection with well
Add the arguments as a list, so pandas doesn't complain about having to set
an index.

.. GENERATED FROM PYTHON SOURCE LINES 68-73

.. code-block:: Python

    weldata = pd.DataFrame()
    weldata["x"] = [0.5]
    weldata["y"] = [0.5]
    weldata["q"] = [qscaled]








.. GENERATED FROM PYTHON SOURCE LINES 74-78

Build
-----

Finally, we build the model.

.. GENERATED FROM PYTHON SOURCE LINES 78-108

.. code-block:: Python


    m = imod.wq.SeawatModel("Henry")
    m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=50.0, bottom=bot, starting_head=1.0)
    m["lpf"] = imod.wq.LayerPropertyFlow(
        k_horizontal=10.0, k_vertical=10.0, specific_storage=0.0
    )
    m["btn"] = imod.wq.BasicTransport(
        icbund=bnd, starting_concentration=35.0, porosity=0.35
    )
    m["adv"] = imod.wq.AdvectionTVD(courant=1.0)
    m["dsp"] = imod.wq.Dispersion(longitudinal=0.1, diffusion_coefficient=1.0e-9)
    m["vdf"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)
    m["wel"] = imod.wq.Well(
        id_name="well", x=weldata["x"], y=weldata["y"], rate=weldata["q"]
    )
    m["pcg"] = imod.wq.PreconditionedConjugateGradientSolver(
        max_iter=150, inner_iter=30, hclose=0.0001, rclose=1.0, relax=0.98, damp=1.0
    )
    m["gcg"] = imod.wq.GeneralizedConjugateGradientSolver(
        max_iter=150,
        inner_iter=30,
        cclose=1.0e-6,
        preconditioner="mic",
        lump_dispersion=True,
    )
    m["oc"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)
    m.create_time_discretization(
        additional_times=pd.date_range("2000-01-01", "2001-01-01", freq="M")
    )








.. GENERATED FROM PYTHON SOURCE LINES 109-110

Now we write the model

.. GENERATED FROM PYTHON SOURCE LINES 110-114

.. code-block:: Python


    modeldir = imod.util.temporary_directory()
    m.write(modeldir, resultdir_is_workdir=True)








.. GENERATED FROM PYTHON SOURCE LINES 115-152

Run
---

You can run the model using the comand prompt and the iMOD-WQ executable.
This is part of the iMOD v5 release, which can be downloaded here:
https://oss.deltares.nl/web/imod/download-imod5 .
This only works on Windows.

To run your model, open up a command prompt
and run the following commands:

.. code-block:: batch

   cd c:\path\to\modeldir
   c:\path\to\imod\folder\iMOD-WQ_V5_3_SVN359_X64R.exe Henry.run

Note that the version name of your executable might differ.

%%
Visualise results
-----------------

After succesfully running the model, you can
plot results as follows:

.. code:: python

   head = imod.idf.open(modeldir / "results/head/*.idf")

   fig, ax = plt.subplots()
   head.plot(yincrease=False, ax=ax)

   conc = imod.idf.open(modeldir / "results/conc/*.idf")

   fig, ax = plt.subplots()
   conc.plot(levels=range(0, 35, 5), yincrease=False, ax=ax)



.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.500 seconds)


.. _sphx_glr_download_examples_imod-wq_Henry-wq.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: Henry-wq.ipynb <Henry-wq.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: Henry-wq.py <Henry-wq.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: Henry-wq.zip <Henry-wq.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_