.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/imod-wq/VerticalInterface.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:Go to the end  to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_imod-wq_VerticalInterface.py: Vertical Interface ================== This 2D examples demonstrates the rotation of an initially vertical interface between fresh and salt water. For a detailed description of this benchmark, see: Bakker, M., Oude Essink, G. H. P., & Langevin, C. D. (2004). The rotating movement of three immiscible fluids - A benchmark problem. Journal of Hydrology, 287 (1-4), 270-278. https://doi.org/10.1016/j.jhydrol.2003.10.007 .. GENERATED FROM PYTHON SOURCE LINES 17-18 We'll start with the usual imports .. GENERATED FROM PYTHON SOURCE LINES 18-27 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import pandas as pd import xarray as xr import imod .. GENERATED FROM PYTHON SOURCE LINES 29-35 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 35-52 .. code-block:: Python nrow = 1 ncol = 80 nlay = 40 dz = 1.0 # 0.0125 dx = 1.0 # 0.0125 dy = -dx # Defining tops and bottoms top1D = xr.DataArray( np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, nlay + 1)}, ("layer") ) bot = top1D - dz top = nlay * dz .. GENERATED FROM PYTHON SOURCE LINES 53-54 Set up ibound, which sets where active cells are (ibound = 1.0) .. GENERATED FROM PYTHON SOURCE LINES 54-69 .. 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"), ) fig, ax = plt.subplots() bnd.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_VerticalInterface_001.png :alt: y = 0.5, dx = 1.0, dy = -1.0 :srcset: /examples/imod-wq/images/sphx_glr_VerticalInterface_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 70-74 Define the icbund, which sets which cells in the solute transport model are active, inactive or constant. We just go for active cells everywhere here .. GENERATED FROM PYTHON SOURCE LINES 74-76 .. code-block:: Python icbund = xr.full_like(bnd, 1) .. GENERATED FROM PYTHON SOURCE LINES 77-82 Boundary Conditions ------------------- Set the constant heads by specifying a negative value in iboud, that is: bnd[index] = -1 .. GENERATED FROM PYTHON SOURCE LINES 82-89 .. code-block:: Python bnd[31, :, 0] = -1 fig, ax = plt.subplots() bnd.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_VerticalInterface_002.png :alt: y = 0.5, dx = 1.0, dy = -1.0 :srcset: /examples/imod-wq/images/sphx_glr_VerticalInterface_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 90-91 Define WEL data .. GENERATED FROM PYTHON SOURCE LINES 91-96 .. code-block:: Python weldata = pd.DataFrame() weldata["x"] = np.full(1, 0.5 * dx) weldata["y"] = np.full(1, 0.5) weldata["q"] = 0.28512 # positive, so it's an injection well .. GENERATED FROM PYTHON SOURCE LINES 97-101 Initial Conditions ------------------ Define the starting concentration .. GENERATED FROM PYTHON SOURCE LINES 101-114 .. code-block:: Python sconc = xr.DataArray( data=np.full((nlay, nrow, ncol), 0.0), coords={ "y": [0.5], "x": np.arange(0.5 * dx, dx * ncol, dx), "layer": np.arange(1, nlay + 1), }, dims=("layer", "y", "x"), ) sconc[:, :, 41:80] = 35.0 .. GENERATED FROM PYTHON SOURCE LINES 115-119 Build ----- Finally, we build the model. .. GENERATED FROM PYTHON SOURCE LINES 119-151 .. code-block:: Python fig, ax = plt.subplots() sconc.plot(y="layer", yincrease=False, ax=ax) # Finally, we build the model m = imod.wq.SeawatModel("VerticalInterface") m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=top, bottom=bot, starting_head=0.0) m["lpf"] = imod.wq.LayerPropertyFlow( k_horizontal=86.4, k_vertical=86.4, specific_storage=0.0 ) m["btn"] = imod.wq.BasicTransport( icbund=icbund, starting_concentration=sconc, porosity=0.1 ) m["adv"] = imod.wq.AdvectionTVD(courant=1.0) m["dsp"] = imod.wq.Dispersion(longitudinal=0.0, diffusion_coefficient=0.0) m["vdf"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71) m["wel"] = imod.wq.Well( id_name="wel", x=weldata["x"], y=weldata["y"], rate=weldata["q"] ) m["pcg"] = imod.wq.PreconditionedConjugateGradientSolver( max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, 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=["2000-01-01", "2000-01-02"]) .. image-sg:: /examples/imod-wq/images/sphx_glr_VerticalInterface_003.png :alt: y = 0.5 :srcset: /examples/imod-wq/images/sphx_glr_VerticalInterface_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 152-153 Now we write the model, including runfile: .. GENERATED FROM PYTHON SOURCE LINES 153-156 .. code-block:: Python modeldir = imod.util.temporary_directory() m.write(modeldir, resultdir_is_workdir=True) .. GENERATED FROM PYTHON SOURCE LINES 157-194 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 VerticalInterface.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.554 seconds) .. _sphx_glr_download_examples_imod-wq_VerticalInterface.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:Download Jupyter notebook: VerticalInterface.ipynb  .. container:: sphx-glr-download sphx-glr-download-python :download:Download Python source code: VerticalInterface.py  .. only:: html .. rst-class:: sphx-glr-signature Gallery generated by Sphinx-Gallery `_