.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/imod-wq/Hydrocoin.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:here  to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_imod-wq_Hydrocoin.py: Hydrocoin ========= A 2D case from the Hydrological Code Intercomparison (Hydrocoin). For more information see: Konikow, L. F., Sanford, W. E., & Campbell, P. J. (1997). Constant-concentration boundary condition: Lessons from the HYDROCOIN variable-density groundwater benchmark problem. Water Resources Research, 33 (10), 2253-2261. https://doi.org/10.1029/97WR01926 .. GENERATED FROM PYTHON SOURCE LINES 14-17 .. code-block:: default import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 18-19 We'll start with the usual imports .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. code-block:: default import numpy as np import pandas as pd import xarray as xr import imod .. GENERATED FROM PYTHON SOURCE LINES 28-34 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 34-43 .. code-block:: default nrow = 1 # number of rows ncol = 45 # number of columns nlay = 76 # number of layers dz = 4.0 dx = 20.0 dy = -dx .. GENERATED FROM PYTHON SOURCE LINES 44-45 Set up tops and bottoms .. GENERATED FROM PYTHON SOURCE LINES 45-51 .. code-block:: default top1D = xr.DataArray( np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, nlay + 1)}, ("layer") ) bot = top1D - dz .. GENERATED FROM PYTHON SOURCE LINES 52-53 Set up ibound, which sets where active cells are (ibound = 1.0). .. GENERATED FROM PYTHON SOURCE LINES 53-66 .. code-block:: default 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"), ) .. GENERATED FROM PYTHON SOURCE LINES 67-68 Set inactive cells by specifying bnd[index] = 0.0 .. GENERATED FROM PYTHON SOURCE LINES 68-75 .. code-block:: default bnd[75, :, 0:15] = 0.0 bnd[75, :, 30:45] = 0.0 fig, ax = plt.subplots() bnd.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_Hydrocoin_001.png :alt: y = 0.5, dx = 20.0, dy = -20.0 :srcset: /examples/imod-wq/images/sphx_glr_Hydrocoin_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 76-81 Boundary Conditions ------------------- Set the constant heads by specifying a negative value in iboud, that is: bnd[index] = -1 .. GENERATED FROM PYTHON SOURCE LINES 81-89 .. code-block:: default bnd[0, :, :] = -1 fig, ax = plt.subplots() bnd.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_Hydrocoin_002.png :alt: y = 0.5, dx = 20.0, dy = -20.0 :srcset: /examples/imod-wq/images/sphx_glr_Hydrocoin_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, need to define the x, y, and pumping rate (q) .. GENERATED FROM PYTHON SOURCE LINES 91-96 .. code-block:: default 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-102 Define the icbund, which sets which cells in the solute transport model are active, inactive or constant. In this case the central 15 cells on the top row have a constant concentration, And, on both sides, the outer 15 cells of the top row are inactive in the transport model. .. GENERATED FROM PYTHON SOURCE LINES 102-122 .. code-block:: default icbund = 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, nlay + 1), "dx": dx, "dy": dy, }, dims=("layer", "y", "x"), ) icbund[75, :, 0:15] = 0.0 icbund[75, :, 30:45] = 0.0 icbund[75, :, 15:30] = -1.0 fig, ax = plt.subplots() icbund.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_Hydrocoin_003.png :alt: y = 0.5, dx = 20.0, dy = -20.0 :srcset: /examples/imod-wq/images/sphx_glr_Hydrocoin_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 123-127 Initial conditions ------------------ Define the starting concentrations .. GENERATED FROM PYTHON SOURCE LINES 127-146 .. code-block:: default 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), "dx": dx, "dy": dy, }, dims=("layer", "y", "x"), ) sconc[75, :, 15:30] = 280.0 fig, ax = plt.subplots() sconc.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_Hydrocoin_004.png :alt: y = 0.5, dx = 20.0, dy = -20.0 :srcset: /examples/imod-wq/images/sphx_glr_Hydrocoin_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 147-148 Define starting heads, these will be inserted in the Basic Flow (BAS) package .. GENERATED FROM PYTHON SOURCE LINES 148-214 .. code-block:: default shd = 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), "dx": dx, "dy": dy, }, dims=("layer", "y", "x"), ) shd[0, :, :] = np.array( [ 10, 9.772727273, 9.545454545, 9.318181818, 9.090909091, 8.863636364, 8.636363636, 8.409090909, 8.181818182, 7.954545455, 7.727272727, 7.5, 7.272727273, 7.045454545, 6.818181818, 6.590909091, 6.363636364, 6.136363636, 5.909090909, 5.681818182, 5.454545455, 5.227272727, 5, 4.772727273, 4.545454545, 4.318181818, 4.090909091, 3.863636364, 3.636363636, 3.409090909, 3.181818182, 2.954545455, 2.727272727, 2.5, 2.272727273, 2.045454545, 1.818181818, 1.590909091, 1.363636364, 1.136363636, 0.909090909, 0.681818182, 0.454545455, 0.227272727, 0.00, ] ) fig, ax = plt.subplots() shd.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_Hydrocoin_005.png :alt: y = 0.5, dx = 20.0, dy = -20.0 :srcset: /examples/imod-wq/images/sphx_glr_Hydrocoin_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 215-219 Hydrogeology ------------ Define horizontal hydraulic conductivity .. GENERATED FROM PYTHON SOURCE LINES 219-237 .. code-block:: default khv = xr.DataArray( data=np.full((nlay, nrow, ncol), 0.847584), coords={ "y": [0.5], "x": np.arange(0.5 * dx, dx * ncol, dx), "layer": np.arange(1, nlay + 1), "dx": dx, "dy": dy, }, dims=("layer", "y", "x"), ) khv[75, :, 15:30] = 0.0008475 fig, ax = plt.subplots() khv.plot(y="layer", yincrease=False, ax=ax) .. image-sg:: /examples/imod-wq/images/sphx_glr_Hydrocoin_006.png :alt: y = 0.5, dx = 20.0, dy = -20.0 :srcset: /examples/imod-wq/images/sphx_glr_Hydrocoin_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 238-242 Build ----- Finally, we build the model. .. GENERATED FROM PYTHON SOURCE LINES 242-271 .. code-block:: default m = imod.wq.SeawatModel("Hydrocoin") m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=304.0, bottom=bot, starting_head=shd) m["lpf"] = imod.wq.LayerPropertyFlow( k_horizontal=khv, k_vertical=khv, specific_storage=0.0 ) m["btn"] = imod.wq.BasicTransport( icbund=icbund, starting_concentration=sconc, porosity=0.2 ) m["adv"] = imod.wq.AdvectionTVD(courant=1.0) m["dsp"] = imod.wq.Dispersion(longitudinal=20.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-01T00:00", "2010-01-01T00:00"]) .. GENERATED FROM PYTHON SOURCE LINES 272-273 Now we write the model, including runfile: .. GENERATED FROM PYTHON SOURCE LINES 273-276 .. code-block:: default modeldir = imod.util.temporary_directory() m.write(modeldir, resultdir_is_workdir=True) .. GENERATED FROM PYTHON SOURCE LINES 277-314 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 Hydrocoin.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 2.142 seconds) .. _sphx_glr_download_examples_imod-wq_Hydrocoin.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:Download Python source code: Hydrocoin.py  .. container:: sphx-glr-download sphx-glr-download-jupyter :download:Download Jupyter notebook: Hydrocoin.ipynb  .. only:: html .. rst-class:: sphx-glr-signature Gallery generated by Sphinx-Gallery _