.. 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 27-29 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 29-35 .. code-block:: Python import numpy as np import xarray as xr import imod .. GENERATED FROM PYTHON SOURCE LINES 36-41 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 41-60 .. 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 61-71 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 71-107 .. 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 108-113 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 113-181 .. 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", csv_output=False, no_ptc=True, 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 182-183 We'll create a new directory in which we will write and run the model. .. GENERATED FROM PYTHON SOURCE LINES 183-187 .. code-block:: Python modeldir = imod.util.temporary_directory() simulation.write(modeldir) .. GENERATED FROM PYTHON SOURCE LINES 188-196 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 196-199 .. code-block:: Python simulation.run() .. GENERATED FROM PYTHON SOURCE LINES 200-204 Open the results ---------------- We'll open the heads (.hds) file. .. GENERATED FROM PYTHON SOURCE LINES 204-210 .. 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 211-213 Visualize the results --------------------- .. GENERATED FROM PYTHON SOURCE LINES 213-217 .. 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 218-226 .. _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.323 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 ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_