.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\imodflow\ImodflowProjectfile.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_imodflow_ImodflowProjectfile.py: Model Creation ============== In this example, we'll create an iMODFLOW model from scratch with complex boundary conditions and horizontal barriers. There are two surface water systems: the outer two edges of the grid feature ditches with a rising stage, whereas the central ditch has a periodic boundary conditions with a summer and winter stage. The model will be written as a projectfile with a set of IDFs containing all the grid information and a .tim file containing the time discretization. .. GENERATED FROM PYTHON SOURCE LINES 20-23 We'll start with the usual imports, supplied with geopandas and shapely to specify vector data for the `hfb` package. .. GENERATED FROM PYTHON SOURCE LINES 23-33 .. code-block:: Python import geopandas as gpd import numpy as np import pandas as pd import xarray as xr from shapely.geometry import LineString import imod import imod.flow as flow .. GENERATED FROM PYTHON SOURCE LINES 34-41 Discretization -------------- We'll start off by creating a model discretization. The model consists of a 3 by 9 by 9 three-dimensional grid. We'll specify the grid first. .. GENERATED FROM PYTHON SOURCE LINES 41-52 .. code-block:: Python shape = nlay, nrow, ncol = 3, 9, 9 dx = 100.0 dy = -100.0 dz = np.array([5, 30, 100]) xmin = 0.0 xmax = dx * ncol ymin = 0.0 ymax = abs(dy) * nrow dims = ("layer", "y", "x") .. GENERATED FROM PYTHON SOURCE LINES 53-55 Next, we'll create the coordinates which set the grid dimensions. .. GENERATED FROM PYTHON SOURCE LINES 55-61 .. code-block:: Python layer = np.arange(1, nlay + 1) y = np.arange(ymax, ymin, dy) + 0.5 * dy x = np.arange(xmin, xmax, dx) + 0.5 * dx coords = {"layer": layer, "y": y, "x": x} .. GENERATED FROM PYTHON SOURCE LINES 62-63 The vertical grid discretization (tops and bottoms) are set with a 1D DataArray. .. GENERATED FROM PYTHON SOURCE LINES 63-70 .. code-block:: Python surface = 0.0 interfaces = np.insert((surface - np.cumsum(dz)), 0, surface) bottom = xr.DataArray(interfaces[1:], coords={"layer": layer}, dims="layer") top = xr.DataArray(interfaces[:-1], coords={"layer": layer}, dims="layer") .. GENERATED FROM PYTHON SOURCE LINES 71-73 We'll have to create a time discretization as well. Create 1 year of monthly timesteps .. GENERATED FROM PYTHON SOURCE LINES 73-75 .. code-block:: Python times = pd.date_range(start="1/1/2018", end="12/1/2018", freq="MS") .. GENERATED FROM PYTHON SOURCE LINES 76-78 We'll create our first 3 dimensional grid here, the `ibound` grid, which sets where active cells are `(ibound = 1.0)` .. GENERATED FROM PYTHON SOURCE LINES 78-81 .. code-block:: Python ibound = xr.DataArray(np.ones(shape), coords=coords, dims=dims) ibound .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 82-87 Hydrogeology ------------ We'll create a very simple hydrogeology, by specifying kh, kva, and sto as constants .. GENERATED FROM PYTHON SOURCE LINES 87-92 .. code-block:: Python kh = 10.0 kva = 1.0 sto = 0.001 .. GENERATED FROM PYTHON SOURCE LINES 93-100 Initial conditions ------------------ We do not put much effort in the creation of the initial conditions in this example, instead we copy the ibounds. This is a 3D grid filled with the value 1, and we can use it as a inital condition as well. .. GENERATED FROM PYTHON SOURCE LINES 100-102 .. code-block:: Python starting_head = ibound.copy() .. GENERATED FROM PYTHON SOURCE LINES 103-117 Boundary conditions ------------------- We will put some more effort in creating some complex boundary conditions. We'll create both two outer ditches with a rising stage, as well as a central ditch with periodic (summer-winter) stage. Rising outer ditches ~~~~~~~~~~~~~~~~~~~~ First, we'll create rising trend, by creating a 1D array ones with the same size as the time dimension, and computing the cumulative sum over it. .. GENERATED FROM PYTHON SOURCE LINES 117-123 .. code-block:: Python trend = np.ones(times[:-1].shape) trend = np.cumsum(trend) trend_da = xr.DataArray(trend, coords={"time": times[:-1]}, dims=["time"]) trend_da .. raw:: html

<xarray.DataArray (time: 11)> Size: 88B
    array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11.])
    Coordinates:
      * time     (time) datetime64[ns] 88B 2018-01-01 2018-02-01 ... 2018-11-01

.. GENERATED FROM PYTHON SOURCE LINES 124-125 Next, we assign values only to edges of model x domain. .. GENERATED FROM PYTHON SOURCE LINES 125-130 .. code-block:: Python is_x_edge = starting_head.x.isin([x[0], x[-1]]) head_edge = starting_head.where(is_x_edge) head_edge .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 131-137 Now let's multiply our 1D DataArray with dimension ``time``, with the static 3D grid with dimension ``layer, y, x``, which xarray automatically broadcasts to a 4D array, with dimensions ``time, layer, y, x`` This finishes our .. GENERATED FROM PYTHON SOURCE LINES 137-141 .. code-block:: Python head_edge_rising = trend_da * head_edge head_edge_rising .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 142-149 Periodic central ditch ~~~~~~~~~~~~~~~~~~~~~~ We'll take only the central column of the grid with (``where``), the rest will be set to ``np.nan``, and from this we'll select only the upper layer, as the ditch will be located only in the upper layer. .. GENERATED FROM PYTHON SOURCE LINES 149-153 .. code-block:: Python is_x_central = starting_head.x == x[4] head_central = starting_head.where(is_x_central).sel(layer=1) .. GENERATED FROM PYTHON SOURCE LINES 154-162 Create period times, we let these times start before the model starts. This is necessary because `iMODFLOW` only forward fills periods. Otherwise, in this case there wouldn't be a periodic boundary condition until april. We will do this by selecting the months april and october, and then subtracting a year .. GENERATED FROM PYTHON SOURCE LINES 162-164 .. code-block:: Python period_times = times[[3, 9]] - np.timedelta64(365, "D") .. GENERATED FROM PYTHON SOURCE LINES 165-166 We are creating a summer and winter level. .. GENERATED FROM PYTHON SOURCE LINES 166-171 .. code-block:: Python periods_da = xr.DataArray([4, 10], coords={"time": period_times}, dims=["time"]) head_periodic = periods_da * head_central head_periodic .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 172-174 Create dictionary to tell `iMOD` which period name corresponds to which date. .. GENERATED FROM PYTHON SOURCE LINES 174-179 .. code-block:: Python timemap = { period_times[0]: "summer", period_times[1]: "winter", } .. GENERATED FROM PYTHON SOURCE LINES 180-187 Wells ~~~~~ Wells are specified as a pandas dataframe. We create a diagonal line of wells through the domain. Because we can. .. GENERATED FROM PYTHON SOURCE LINES 187-198 .. code-block:: Python wel_df = pd.DataFrame() wel_df["id_name"] = np.arange(len(x)).astype(str) wel_df["x"] = x wel_df["y"] = y wel_df["rate"] = dx * dy * -1 * 0.5 wel_df["time"] = np.tile(times[:-1], 2)[: len(x)] wel_df["layer"] = 2 wel_df .. raw:: html

	id_name	x	y	rate	time	layer
0	0	50.0	850.0	5000.0	2018-01-01	2
1	1	150.0	750.0	5000.0	2018-02-01	2
2	2	250.0	650.0	5000.0	2018-03-01	2
3	3	350.0	550.0	5000.0	2018-04-01	2
4	4	450.0	450.0	5000.0	2018-05-01	2
5	5	550.0	350.0	5000.0	2018-06-01	2
6	6	650.0	250.0	5000.0	2018-07-01	2
7	7	750.0	150.0	5000.0	2018-08-01	2
8	8	850.0	50.0	5000.0	2018-09-01	2

.. GENERATED FROM PYTHON SOURCE LINES 199-203 Horizontal Flow Barrier ----------------------- Create barriers between ditches in layer 1 and 2 (but not 3). .. GENERATED FROM PYTHON SOURCE LINES 203-207 .. code-block:: Python line1 = LineString([(x[2], ymin), (x[2], ymax)]) line2 = LineString([(x[7], ymin), (x[7], ymax)]) .. GENERATED FROM PYTHON SOURCE LINES 208-209 We'll have to repeat each line for each layer .. GENERATED FROM PYTHON SOURCE LINES 209-212 .. code-block:: Python lines = np.array([line1, line2, line1, line2], dtype="object") hfb_layers = np.array([1, 1, 2, 2]) .. GENERATED FROM PYTHON SOURCE LINES 213-214 We can specify names for our own bookkeeping .. GENERATED FROM PYTHON SOURCE LINES 214-225 .. code-block:: Python id_name = ["left_upper", "right_upper", "left_lower", "right_lower"] # The hfb has to specified as a geopandas `GeoDataFrame # `_ hfb_gdf = gpd.GeoDataFrame( geometry=lines, data={"id_name": id_name, "layer": hfb_layers, "resistance": 100.0} ) hfb_gdf .. raw:: html

	id_name	layer	resistance	geometry
0	left_upper	1	100.0	LINESTRING (250 0, 250 900)
1	right_upper	1	100.0	LINESTRING (750 0, 750 900)
2	left_lower	2	100.0	LINESTRING (250 0, 250 900)
3	right_lower	2	100.0	LINESTRING (750 0, 750 900)

.. GENERATED FROM PYTHON SOURCE LINES 226-230 Build ----- Finally, we build the model. .. GENERATED FROM PYTHON SOURCE LINES 230-246 .. code-block:: Python m = flow.ImodflowModel("my_first_imodflow_model") m["pcg"] = flow.PreconditionedConjugateGradientSolver() m["bnd"] = flow.Boundary(ibound) m["top"] = flow.Top(top) m["bottom"] = flow.Bottom(bottom) m["khv"] = flow.HorizontalHydraulicConductivity(kh) m["kva"] = flow.VerticalAnisotropy(kva) m["sto"] = flow.StorageCoefficient(sto) m["shd"] = flow.StartingHead(starting_head) m["chd"] = flow.ConstantHead(head=head_edge_rising) .. GENERATED FROM PYTHON SOURCE LINES 247-249 Create periodic boundary condition and specify it as a periodic stress package. .. GENERATED FROM PYTHON SOURCE LINES 249-252 .. code-block:: Python m["ghb"] = flow.GeneralHeadBoundary(head=head_periodic, conductance=10.0) m["ghb"].periodic_stress(timemap) .. GENERATED FROM PYTHON SOURCE LINES 253-254 We can specify a second stress package as follows: .. GENERATED FROM PYTHON SOURCE LINES 254-262 .. code-block:: Python m["ghb2"] = flow.GeneralHeadBoundary(head=head_periodic + 10.0, conductance=1.0) # You also need to specify periodic stresses for second system. m["ghb2"].periodic_stress(timemap) m["wel"] = flow.Well(**wel_df) m["hfb"] = flow.HorizontalFlowBarrier(**hfb_gdf) .. GENERATED FROM PYTHON SOURCE LINES 263-264 Specify output control; -1 specifies to save the output of interest for all layers. .. GENERATED FROM PYTHON SOURCE LINES 264-267 .. code-block:: Python m["oc"] = flow.OutputControl(save_head=-1, save_flux=-1) .. GENERATED FROM PYTHON SOURCE LINES 268-269 imod-python wants you to provide the model endtime to your time_discretization! .. GENERATED FROM PYTHON SOURCE LINES 269-271 .. code-block:: Python m.create_time_discretization(additional_times=times[-1]) .. GENERATED FROM PYTHON SOURCE LINES 272-275 Now we write the model Writes both .IDFs as well as projectfile, an inifile, and a .tim file that contains the time discretization. .. GENERATED FROM PYTHON SOURCE LINES 275-279 .. code-block:: Python modeldir = imod.util.temporary_directory() m.write(directory=modeldir) .. GENERATED FROM PYTHON SOURCE LINES 280-302 Run --- You can run the model using the comand prompt and the iMOD executables. This is part of the iMOD v5 release, which can be downloaded here: https://oss.deltares.nl/web/imod/download-imod5 . iMOD 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_v5_3_X64R.EXE config_run.ini c:\path\to\imod\folder\iMODFLOW_V5_3_METASWAP_SVN1977_X64R.EXE my_first_imodflow_model.nam .. note:: iMODFLOW requires you to change directory into the model directory before calling its executable. Otherwise the program will throw an error. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.733 seconds) .. _sphx_glr_download_examples_imodflow_ImodflowProjectfile.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ImodflowProjectfile.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ImodflowProjectfile.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ImodflowProjectfile.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_