.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\mf6\example_models.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_example_models.py: Example models ============== This source file contains functions that create a simulation that can be used in examples that are not focused on building a simulation, but on doing something with it (such as regridding). .. GENERATED FROM PYTHON SOURCE LINES 10-296 .. code-block:: Python import numpy as np import xarray as xr import xugrid as xu import imod def create_twri_simulation() -> imod.mf6.Modflow6Simulation: """There is a separate example contained in :doc:`/examples/mf6/ex01_twri` that you should look at if you are interested in the model building. The TWRI model has 3 layers and contains wells, a drain and recharge. Geometrically it is rectangular with prescribed head on some of the boundaries. Conductivity is highly anisotropic but constant in each layer. Simulation is steady state. """ 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 = 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, "dx": dx, "dy": -dy} # %% # 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. 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 screen_layer = [2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] # we set the screen top and bottoms such that each well falls in one layer and is long enough not to be filtered out perforation_length = 50 delta_z = 0.1 screen_bottom = bottom[screen_layer] + delta_z screen_top = screen_bottom + delta_z + perforation_length # we compute the x and y cooordinates of the wells based on the row and column indices presented in the original twri model well_y = ( ymax - np.array( [ 5.0, 4.0, 6.0, 9.0, 9.0, 9.0, 9.0, 11.0, 11.0, 11.0, 11.0, 13.0, 13.0, 13.0, 13.0, ] ) * abs(dy) + dy / 2 ) well_x = ( np.array( [ 11.0, 6.0, 12.0, 8.0, 10.0, 12.0, 14.0, 8.0, 10.0, 12.0, 14.0, 8.0, 10.0, 12.0, 14.0, ] ) * dx - dx / 2 ) well_rate = [-5.0] * 15 # 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",)) # %% # Build 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. 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=well_x, y=well_y, screen_top=screen_top, screen_bottom=screen_bottom, rate=well_rate, 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"] ) return simulation def create_circle_simulation(): """ There is a separate example contained in :doc:`/examples/mf6/circle` that you should look at if you are interested in the model building. The circle model uses an unstructured grid. It has 2 layers of constant thickness. In conductivity, it is isotropic, and constant in space. Boundary conditions include recharge and constant head. It is a transient simulation. """ # %% # Create a mesh # ------------- grid = imod.data.circle() nface = grid.n_face nlayer = 2 idomain = xu.UgridDataArray( xr.DataArray( np.ones((nlayer, nface), dtype=np.int32), coords={"layer": [1, 2]}, dims=["layer", grid.face_dimension], ), grid=grid, ) # Create model and packages icelltype = xu.full_like(idomain, 0) k = xu.full_like(idomain, 1.0, dtype=float) k33 = k.copy() rch_rate = xu.full_like(idomain.sel(layer=1), 0.001, dtype=float) bottom = idomain * xr.DataArray([5.0, 0.0], dims=["layer"]) chd_location = xu.zeros_like( idomain.sel(layer=2), dtype=bool ).ugrid.binary_dilation(border_value=True) constant_head = xu.full_like(idomain.sel(layer=2), 1.0, dtype=float).where( chd_location ) gwf_model = imod.mf6.GroundwaterFlowModel() gwf_model["disv"] = imod.mf6.VerticesDiscretization( top=10.0, bottom=bottom, idomain=idomain ) gwf_model["chd"] = imod.mf6.ConstantHead( constant_head, 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, save_flows=True, save_specific_discharge=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) # Create simulation simulation = imod.mf6.Modflow6Simulation("circle") simulation["GWF_1"] = gwf_model 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, ) simulation.create_time_discretization(additional_times=["2000-01-01", "2000-01-02"]) return simulation .. _sphx_glr_download_examples_mf6_example_models.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_models.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_models.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_models.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_