.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\mf6\transport_2d.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_transport_2d.py: Transport 2d example ==================== The simulation shown here comes from the 1999 MT3DMS report, p 138: Two-Dimensional Transport in a Uniform Flow of solute injected continuously from a point source in a steady-state uniform flow field. In this example, we build up the model, and the we run the model as is. Next, we split the model in 4 partitions and run that as well. Finally, we show that the difference in outcome for the partitioned and unpartitioned models is small. MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User's Guide .. GENERATED FROM PYTHON SOURCE LINES 20-30 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import pandas as pd import xarray as xr import imod from imod.mf6.multimodel.partition_generator import get_label_array from imod.typing.grid import nan_like, zeros_like .. GENERATED FROM PYTHON SOURCE LINES 31-32 Set some grid dimensions .. GENERATED FROM PYTHON SOURCE LINES 32-42 .. code-block:: Python nlay = 1 # Number of layers nrow = 31 # Number of rows ncol = 46 # Number of columns delr = 10.0 # Column width ($m$) delc = 10.0 # Row width ($m$) delz = 10.0 # Layer thickness ($m$) shape = (nlay, nrow, ncol) top = 10.0 dims = ("layer", "y", "x") .. GENERATED FROM PYTHON SOURCE LINES 43-44 Construct the "idomain" array, and then the discretization package which represents the model grid. .. GENERATED FROM PYTHON SOURCE LINES 44-56 .. code-block:: Python y = np.arange(delr * nrow, 0, -delr) x = np.arange(0, delc * ncol, delc) coords = {"layer": [1], "y": y, "x": x, "dx": delc, "dy": -delr} idomain = xr.DataArray(np.ones(shape, dtype=int), coords=coords, dims=dims) bottom = xr.DataArray([0.0], {"layer": [1]}, ("layer",)) gwf_model = imod.mf6.GroundwaterFlowModel() gwf_model["dis"] = imod.mf6.StructuredDiscretization( top=10.0, bottom=bottom, idomain=idomain ) .. GENERATED FROM PYTHON SOURCE LINES 57-60 Construct the other flow packages. Flow is steady-state in this simulation, meaning specific storage is set to zero. .. GENERATED FROM PYTHON SOURCE LINES 60-72 .. code-block:: Python gwf_model["sto"] = imod.mf6.SpecificStorage( specific_storage=0.0, specific_yield=0.0, transient=False, convertible=0, ) gwf_model["npf"] = imod.mf6.NodePropertyFlow( icelltype=zeros_like(idomain), k=1.0, save_flows=True, save_specific_discharge=True ) gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all") gwf_model["ic"] = imod.mf6.InitialConditions(start=10.0) .. GENERATED FROM PYTHON SOURCE LINES 73-76 Set up the boundary conditions. We have: 2 constant head boundaries at the left and right, chosen so that the velocity is 1/3 m/day and a well that injects 1 m3 per day, with a concentration of 1000 .. GENERATED FROM PYTHON SOURCE LINES 76-116 .. code-block:: Python Lx = 460 v = 1.0 / 3.0 prsity = 0.3 q = v * prsity h1 = q * Lx chd_field = nan_like(idomain) chd_field.values[0, :, 0] = h1 chd_field.values[0, :, -1] = 0.0 chd_concentration = nan_like(idomain) chd_concentration.values[0, :, 0] = 0.0 chd_concentration.values[0, :, -1] = 0.0 chd_concentration = chd_concentration.expand_dims(species=["Au"]) gwf_model["chd"] = imod.mf6.ConstantHead( chd_field, concentration=chd_concentration, print_input=True, print_flows=True, save_flows=True, ) injection_concentration = xr.DataArray( [[1000.0]], coords={ "species": ["Au"], "index": [0], }, dims=("species", "index"), ) gwf_model["wel"] = imod.mf6.Well( x=[150.0], y=[150.0], screen_top=[10.0], screen_bottom=[0.0], rate=[1.0], concentration=injection_concentration, concentration_boundary_type="aux", ) .. GENERATED FROM PYTHON SOURCE LINES 117-121 Now construct the transport simulation. The flow boundaries already have inflow concentration data associated, so the transport boundaries can be imported using the ssm package, and the rest of the transport model definition is straightforward. .. GENERATED FROM PYTHON SOURCE LINES 121-141 .. code-block:: Python tpt_model = imod.mf6.GroundwaterTransportModel() tpt_model["ssm"] = imod.mf6.SourceSinkMixing.from_flow_model( gwf_model, species="Au", save_flows=True ) tpt_model["adv"] = imod.mf6.AdvectionUpstream() tpt_model["dsp"] = imod.mf6.Dispersion( diffusion_coefficient=0.0, longitudinal_horizontal=10.0, transversal_horizontal1=3.0, xt3d_off=False, xt3d_rhs=False, ) tpt_model["mst"] = imod.mf6.MobileStorageTransfer( porosity=0.3, ) tpt_model["ic"] = imod.mf6.InitialConditions(start=0.0) tpt_model["oc"] = imod.mf6.OutputControl(save_concentration="all", save_budget="last") tpt_model["dis"] = gwf_model["dis"] .. GENERATED FROM PYTHON SOURCE LINES 142-145 Create a simulation and add the flow and transport models to it. Then define some ims packages: 1 for every type of model. Finally create 365 time steps of 1 day each. .. GENERATED FROM PYTHON SOURCE LINES 145-187 .. code-block:: Python simulation = imod.mf6.Modflow6Simulation("ex01-twri") simulation["GWF_1"] = gwf_model simulation["TPT_1"] = tpt_model simulation["flow_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["transport_solver"] = imod.mf6.Solution( modelnames=["TPT_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="bicgstab", scaling_method=None, reordering_method=None, relaxation_factor=0.97, ) # Collect time discretization duration = pd.to_timedelta("365d") start = pd.to_datetime("2002-01-01") simulation.create_time_discretization(additional_times=[start, start + duration]) simulation["time_discretization"]["n_timesteps"] = 365 modeldir = imod.util.temporary_directory() simulation.write(modeldir, binary=False) .. GENERATED FROM PYTHON SOURCE LINES 188-190 To split the model in 4 partitions, we must create a label array. We use the utility function ``get_label_array`` for that. .. GENERATED FROM PYTHON SOURCE LINES 190-196 .. code-block:: Python label_array = get_label_array(simulation, 4) fig, ax = plt.subplots() label_array.plot(ax=ax) split_simulation = simulation.split(label_array) .. image-sg:: /examples/mf6/images/sphx_glr_transport_2d_001.png :alt: layer = 1, dx = 10.0, dy = -10.0 :srcset: /examples/mf6/images/sphx_glr_transport_2d_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 197-198 Run the unsplit model and load the simulation results. .. GENERATED FROM PYTHON SOURCE LINES 198-201 .. code-block:: Python simulation.run() conc = simulation.open_concentration() .. GENERATED FROM PYTHON SOURCE LINES 202-203 Run the split model and load the simulation results. .. GENERATED FROM PYTHON SOURCE LINES 203-210 .. code-block:: Python split_modeldir = modeldir / "split" split_simulation.write(split_modeldir, binary=False) split_simulation.run() split_conc = split_simulation.open_concentration()["concentration"] fig, ax = plt.subplots() split_conc.isel(layer=0, time=364).plot.contourf(ax=ax, levels=[0.1, 1, 10]) .. image-sg:: /examples/mf6/images/sphx_glr_transport_2d_002.png :alt: dx = 10.0, dy = -10.0, layer = 1, time = 365.0 :srcset: /examples/mf6/images/sphx_glr_transport_2d_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 211-213 Compute the difference between the split and unsplit simulation results for transport at the end of the simulation, and print them .. GENERATED FROM PYTHON SOURCE LINES 213-216 .. code-block:: Python diff = abs(conc - split_conc) fig, ax = plt.subplots() diff.isel(layer=0, time=364).plot.contourf(ax=ax) .. image-sg:: /examples/mf6/images/sphx_glr_transport_2d_003.png :alt: dx = 10.0, dy = -10.0, layer = 1, time = 365.0 :srcset: /examples/mf6/images/sphx_glr_transport_2d_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 0.979 seconds) .. _sphx_glr_download_examples_mf6_transport_2d.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: transport_2d.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: transport_2d.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: transport_2d.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_