.. 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 Click :ref:here  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:: default 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:: default 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:: default 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:: default 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:: default 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:: default ibound = xr.DataArray(np.ones(shape), coords=coords, dims=dims) ibound .. raw:: html
<xarray.DataArray (layer: 3, y: 9, x: 9)>
array([[[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.]],

[[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.]],

[[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1.]]])
Coordinates:
* layer    (layer) int64 1 2 3
* y        (y) float64 850.0 750.0 650.0 550.0 450.0 350.0 250.0 150.0 50.0
* x        (x) float64 50.0 150.0 250.0 350.0 450.0 550.0 650.0 750.0 850.0

.. 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:: default 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:: default 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:: default 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)>
array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11.])
Coordinates:
* time     (time) datetime64[ns] 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:: default is_x_edge = starting_head.x.isin([x[0], x[-1]]) head_edge = starting_head.where(is_x_edge) head_edge .. raw:: html
<xarray.DataArray (layer: 3, y: 9, x: 9)>
array([[[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.]],

[[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.]],

[[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.],
[ 1., nan, nan, nan, nan, nan, nan, nan,  1.]]])
Coordinates:
* layer    (layer) int64 1 2 3
* y        (y) float64 850.0 750.0 650.0 550.0 450.0 350.0 250.0 150.0 50.0
* x        (x) float64 50.0 150.0 250.0 350.0 450.0 550.0 650.0 750.0 850.0

.. 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:: default head_edge_rising = trend_da * head_edge head_edge_rising .. raw:: html
<xarray.DataArray (time: 11, layer: 3, y: 9, x: 9)>
array([[[[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
...,
[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.]],

[[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
...,
[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.]],

[[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
[ 1., nan, nan, ..., nan, nan,  1.],
...,
...
...,
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.]],

[[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.],
...,
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.]],

[[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.],
...,
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.],
[11., nan, nan, ..., nan, nan, 11.]]]])
Coordinates:
* time     (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-11-01
* layer    (layer) int64 1 2 3
* y        (y) float64 850.0 750.0 650.0 550.0 450.0 350.0 250.0 150.0 50.0
* x        (x) float64 50.0 150.0 250.0 350.0 450.0 550.0 650.0 750.0 850.0

.. 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:: default 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:: default 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:: default periods_da = xr.DataArray([4, 10], coords={"time": period_times}, dims=["time"]) head_periodic = periods_da * head_central head_periodic .. raw:: html
<xarray.DataArray (time: 2, y: 9, x: 9)>
array([[[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan],
[nan, nan, nan, nan,  4., nan, nan, nan, nan]],

[[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan],
[nan, nan, nan, nan, 10., nan, nan, nan, nan]]])
Coordinates:
* time     (time) datetime64[ns] 2017-04-01 2017-10-01
layer    int64 1
* y        (y) float64 850.0 750.0 650.0 550.0 450.0 350.0 250.0 150.0 50.0
* x        (x) float64 50.0 150.0 250.0 350.0 450.0 550.0 650.0 750.0 850.0

.. 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:: default 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:: default 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:: default 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:: default lines = np.array([line1, line2, line1, line2], dtype="object") hfb_layers = np.array([1, 1, 2, 2]) .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/conda/envs/imod/lib/python3.10/site-packages/sphinx_gallery/gen_rst.py:600: ShapelyDeprecationWarning: The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead. exec(self.code, self.fake_main.__dict__) /builds/deltares/imod/imod-python/examples/imodflow/ImodflowProjectfile.py:209: FutureWarning: The input object of type 'LineString' is an array-like implementing one of the corresponding protocols (__array__, __array_interface__ or __array_struct__); but not a sequence (or 0-D). In the future, this object will be coerced as if it was first converted using np.array(obj). To retain the old behaviour, you have to either modify the type 'LineString', or assign to an empty array created with np.empty(correct_shape, dtype=object). lines = np.array([line1, line2, line1, line2], dtype="object") .. GENERATED FROM PYTHON SOURCE LINES 213-214 We can specify names for our own bookkeeping .. GENERATED FROM PYTHON SOURCE LINES 214-225 .. code-block:: default 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=dict(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.000 0.000, 250.000 900.000)
1 right_upper 1 100.0 LINESTRING (750.000 0.000, 750.000 900.000)
2 left_lower 2 100.0 LINESTRING (250.000 0.000, 250.000 900.000)
3 right_lower 2 100.0 LINESTRING (750.000 0.000, 750.000 900.000)

.. GENERATED FROM PYTHON SOURCE LINES 226-230 Build ----- Finally, we build the model. .. GENERATED FROM PYTHON SOURCE LINES 230-246 .. code-block:: default 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:: default 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:: default 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:: default 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:: default 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:: default 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.342 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-python :download:Download Python source code: ImodflowProjectfile.py  .. container:: sphx-glr-download sphx-glr-download-jupyter :download:Download Jupyter notebook: ImodflowProjectfile.ipynb  .. only:: html .. rst-class:: sphx-glr-signature Gallery generated by Sphinx-Gallery _