.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\mf6\ex01_twri.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_mf6_ex01_twri.py: TWRI ==== This example has been converted from the `MODFLOW6 Example problems`_. See the `description`_ and the `notebook`_ which uses `FloPy`_ to setup the model. This example is a modified version of the original MODFLOW example ("`Techniques of Water-Resources Investigation`_" (TWRI)) described in (`McDonald & Harbaugh, 1988`_) and duplicated in (`Harbaugh & McDonald, 1996`_). This problem is also is distributed with MODFLOW-2005 (`Harbaugh, 2005`_). The problem has been modified from a quasi-3D problem, where confining beds are not explicitly simulated, to an equivalent three-dimensional problem. In overview, we'll set the following steps: * Create a structured grid for a rectangular geometry. * Create the xarray DataArrays containg the MODFLOW6 parameters. * Feed these arrays into the imod mf6 classes. * Write to modflow6 files. * Run the model. * Open the results back into DataArrays. * Visualize the results. .. GENERATED FROM PYTHON SOURCE LINES 28-30 We'll start with the usual imports. As this is an simple (synthetic) structured model, we can make due with few packages. .. GENERATED FROM PYTHON SOURCE LINES 30-36 .. code-block:: Python import numpy as np import xarray as xr import imod .. GENERATED FROM PYTHON SOURCE LINES 37-42 Create grid coordinates ----------------------- The first steps consist of setting up the grid -- first the number of layer, rows, and columns. Cell sizes are constant throughout the model. .. GENERATED FROM PYTHON SOURCE LINES 42-61 .. code-block:: Python nlay = 3 nrow = 15 ncol = 15 shape = (nlay, nrow, ncol) dx = 5000.0 dy = -5000.0 xmin = 0.0 xmax = dx * ncol ymin = 0.0 ymax = abs(dy) * nrow dims = ("layer", "y", "x") layer = np.array([1, 2, 3]) 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-72 Create DataArrays ----------------- Now that we have the grid coordinates setup, we can start defining model parameters. The model is characterized by: * a constant head boundary on the left * a single drain in the center left of the model * uniform recharge on the top layer * a number of wells scattered throughout the model. .. GENERATED FROM PYTHON SOURCE LINES 72-108 .. code-block:: Python idomain = xr.DataArray(np.ones(shape, dtype=int), coords=coords, dims=dims) bottom = xr.DataArray([-200.0, -300.0, -450.0], {"layer": layer}, ("layer",)) # Constant head constant_head = xr.full_like(idomain, np.nan, dtype=float).sel(layer=[1, 2]) constant_head[..., 0] = 0.0 # Drainage elevation = xr.full_like(idomain.sel(layer=1), np.nan, dtype=float) conductance = xr.full_like(idomain.sel(layer=1), np.nan, dtype=float) elevation[7, 1:10] = np.array([0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0]) conductance[7, 1:10] = 1.0 # Recharge rch_rate = xr.full_like(idomain.sel(layer=1), 3.0e-8, dtype=float) # Well # fmt: off wells_x = [52500.0, 27500.0, 57500.0, 37500.0, 47500.0, 57500.0, 67500.0, 37500.0, 47500.0, 57500.0, 67500.0, 37500.0, 47500.0, 57500.0, 67500.0, ] wells_y = [52500.0, 57500.0, 47500.0, 32500.0, 32500.0, 32500.0, 32500.0, 22500.0, 22500.0, 22500.0, 22500.0, 12500.0, 12500.0, 12500.0, 12500.0, ] screen_top = [-300.0, -200.0, -200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, ] screen_bottom = [-450.0, -300.0, -300.0, -200.0, -200.0, -200.0, -200.0, -200.0, -200.0, -200.0, -200.0, -200.0, -200.0, -200.0, -200.0, ] rate_wel = [-5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, ] # fmt: on # Node properties icelltype = xr.DataArray([1, 0, 0], {"layer": layer}, ("layer",)) k = xr.DataArray([1.0e-3, 1.0e-4, 2.0e-4], {"layer": layer}, ("layer",)) k33 = xr.DataArray([2.0e-8, 2.0e-8, 2.0e-8], {"layer": layer}, ("layer",)) .. GENERATED FROM PYTHON SOURCE LINES 109-114 Write the model --------------- The first step is to define an empty model, the parameters and boundary conditions are added in the form of the familiar MODFLOW packages. .. GENERATED FROM PYTHON SOURCE LINES 114-180 .. code-block:: Python gwf_model = imod.mf6.GroundwaterFlowModel() gwf_model["dis"] = imod.mf6.StructuredDiscretization( top=200.0, bottom=bottom, idomain=idomain ) gwf_model["chd"] = imod.mf6.ConstantHead( constant_head, print_input=True, print_flows=True, save_flows=True ) gwf_model["drn"] = imod.mf6.Drainage( elevation=elevation, conductance=conductance, print_input=True, print_flows=True, save_flows=True, ) gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0) gwf_model["npf"] = imod.mf6.NodePropertyFlow( icelltype=icelltype, k=k, k33=k33, variable_vertical_conductance=True, dewatered=True, perched=True, save_flows=True, ) gwf_model["sto"] = imod.mf6.SpecificStorage( specific_storage=1.0e-5, specific_yield=0.15, transient=False, convertible=0, ) gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all") gwf_model["rch"] = imod.mf6.Recharge(rch_rate) gwf_model["wel"] = imod.mf6.Well( x=wells_x, y=wells_y, screen_top=screen_top, screen_bottom=screen_bottom, rate=rate_wel, minimum_k=0.0001, ) # Attach it to a simulation simulation = imod.mf6.Modflow6Simulation("ex01-twri") simulation["GWF_1"] = gwf_model # Define solver settings simulation["solver"] = imod.mf6.Solution( modelnames=["GWF_1"], print_option="summary", outer_dvclose=1.0e-4, outer_maximum=500, under_relaxation=None, inner_dvclose=1.0e-4, inner_rclose=0.001, inner_maximum=100, linear_acceleration="cg", scaling_method=None, reordering_method=None, relaxation_factor=0.97, ) # Collect time discretization simulation.create_time_discretization( additional_times=["2000-01-01", "2000-01-02", "2000-01-03", "2000-01-04"] ) .. GENERATED FROM PYTHON SOURCE LINES 181-182 We'll create a new directory in which we will write and run the model. .. GENERATED FROM PYTHON SOURCE LINES 182-186 .. code-block:: Python modeldir = imod.util.temporary_directory() simulation.write(modeldir) .. GENERATED FROM PYTHON SOURCE LINES 187-195 Run the model ------------- .. note:: The following lines assume the ``mf6`` executable is available on your PATH. :ref:`The Modflow 6 examples introduction ` shortly describes how to add it to yours. .. GENERATED FROM PYTHON SOURCE LINES 195-198 .. code-block:: Python simulation.run() .. GENERATED FROM PYTHON SOURCE LINES 199-203 Open the results ---------------- We'll open the heads (.hds) file. .. GENERATED FROM PYTHON SOURCE LINES 203-209 .. code-block:: Python head = imod.mf6.open_hds( modeldir / "GWF_1/GWF_1.hds", modeldir / "GWF_1/dis.dis.grb", ) .. GENERATED FROM PYTHON SOURCE LINES 210-212 Visualize the results --------------------- .. GENERATED FROM PYTHON SOURCE LINES 212-216 .. code-block:: Python head.isel(layer=0, time=0).plot.contourf() .. image-sg:: /examples/mf6/images/sphx_glr_ex01_twri_001.png :alt: dx = 5e+03, dy = -5e+03, layer = 1, time = 1.0 :srcset: /examples/mf6/images/sphx_glr_ex01_twri_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 217-225 .. _MODFLOW6 example problems: https://github.com/MODFLOW-USGS/modflow6-examples .. _description: https://modflow6-examples.readthedocs.io/en/master/_examples/ex-gwf-twri.html .. _notebook: https://github.com/MODFLOW-USGS/modflow6-examples/tree/master/notebooks/ex-gwf-twri.ipynb .. _Techniques of Water-Resources Investigation: https://pubs.usgs.gov/twri/twri7-c1/ .. _McDonald & Harbaugh, 1988: https://pubs.er.usgs.gov/publication/twri06A1 .. _Harbaugh & McDonald, 1996: https://pubs.er.usgs.gov/publication/ofr96485 .. _Harbaugh, 2005: https://pubs.er.usgs.gov/publication/tm6A16 .. _FloPy: https://github.com/modflowpy/flopy .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.548 seconds) .. _sphx_glr_download_examples_mf6_ex01_twri.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ex01_twri.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ex01_twri.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ex01_twri.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_