.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\mf6\circle_transport.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_circle_transport.py: Freshwater lens (circle) ======================== This example illustrates how to setup a very simple unstructured groundwater transport model using the ``imod`` package and associated packages. In overview, we'll set the following steps: * Setting up the flow model, just like in the circle.py example * set up the transport model * Run the simulation. * Visualize the results. .. GENERATED FROM PYTHON SOURCE LINES 15-17 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 19-20 We'll start with the following imports: .. GENERATED FROM PYTHON SOURCE LINES 20-29 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import pandas as pd import xarray as xr import xugrid as xu import imod .. GENERATED FROM PYTHON SOURCE LINES 30-32 Parameters ---------- .. GENERATED FROM PYTHON SOURCE LINES 32-40 .. code-block:: Python porosity = 0.10 max_concentration = 35.0 min_concentration = 0.0 max_density = 1025.0 min_density = 1000.0 k_value = 10.0 .. GENERATED FROM PYTHON SOURCE LINES 41-49 Create a mesh ------------- The first steps consists of generating a mesh. In this example, we'll use data included with iMOD Python for a circular mesh. Note that this is a `Ugrid2D object. `_ For more information on working with unstructured grids see the `Xugrid documentation `_ .. GENERATED FROM PYTHON SOURCE LINES 49-56 .. code-block:: Python grid_triangles = imod.data.circle() fig, ax = plt.subplots() xu.plot.line(grid_triangles) ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_transport_001.png :alt: circle transport :srcset: /examples/mf6/images/sphx_glr_circle_transport_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 57-66 However a triangular grid has the issue that the direction of the fluxes between cell centres is not perpendicular to the cell vertices. The default formulation of Modflow 6 does not account for this, which causes mass balance errors. The XT3D formulation is able to account for this, but the last version of Modflow 6 (6.3 at time of writing) does not support this in combination with the Buoyancy package and using XT3D comes with an extra, significant computational burden. It is therefore easier to use a voronoi grid, for which `Xugrid `_ has a very convenient method. .. GENERATED FROM PYTHON SOURCE LINES 66-73 .. code-block:: Python grid = grid_triangles.tesselate_centroidal_voronoi() fig, ax = plt.subplots() xu.plot.line(grid) ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_transport_002.png :alt: circle transport :srcset: /examples/mf6/images/sphx_glr_circle_transport_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 74-75 Create arrays .. GENERATED FROM PYTHON SOURCE LINES 75-96 .. code-block:: Python nface = grid.n_face nlayer = 15 layer = np.arange(nlayer, dtype=int) + 1 idomain = xu.UgridDataArray( xr.DataArray( np.ones((nlayer, nface), dtype=np.int32), coords={"layer": layer}, dims=["layer", grid.face_dimension], ), grid=grid, ) icelltype = xu.full_like(idomain, 0) k = xu.full_like(idomain, k_value, dtype=float) k33 = k.copy() top = 0.0 bottom = xr.DataArray(top - (layer * 10.0), dims=["layer"]) .. GENERATED FROM PYTHON SOURCE LINES 97-104 Recharge -------- We need a recharge rate for the fluid and a recharge rate for the solute. The fluid recharge rate is volumetric and per unit area, so the unit is length/time. The solute recharge rate is the concentration of solute in the recharge, and has concentration units. .. GENERATED FROM PYTHON SOURCE LINES 104-110 .. code-block:: Python rch_rate = xu.full_like(idomain.sel(layer=1), 0.001, dtype=float) rch_concentration = xu.full_like(rch_rate, min_concentration) rch_concentration = rch_concentration.expand_dims(species=["salinity"]) .. GENERATED FROM PYTHON SOURCE LINES 111-119 Unlike a recharge boundary, with a prescribed head boundary we don't know a priori whether water will flow in over the boundary or leave across the boundary. If water flows into the model over the boundary, it carries a prescribed solute concentration. If it leaves, it leaves with the concentration that was computed for the cell. In this example we set the prescribed head value to 0.0 and the external concentration to 35.0 as well. The boundary only operates on the top layer. .. GENERATED FROM PYTHON SOURCE LINES 119-137 .. code-block:: Python chd_location = xu.zeros_like(idomain.sel(layer=1), dtype=bool).ugrid.binary_dilation( border_value=True ) constant_head = xu.full_like(idomain, 0.0, dtype=float).where(chd_location) # Approximate face area face_area = (1000.0 / 6) ** 2 * 0.5 conductance = xu.full_like(idomain, face_area * k_value, dtype=float).where( chd_location ) constant_concentration = xu.full_like(constant_head, max_concentration).where( chd_location ) constant_concentration = constant_concentration.expand_dims(species=["salinity"]) .. GENERATED FROM PYTHON SOURCE LINES 138-142 Add flow model to simulation ---------------------------- See the circle.py example for more information. .. GENERATED FROM PYTHON SOURCE LINES 142-190 .. code-block:: Python gwf_model = imod.mf6.GroundwaterFlowModel() gwf_model["disv"] = imod.mf6.VerticesDiscretization( top=top, bottom=bottom, idomain=idomain ) gwf_model["ghb"] = imod.mf6.GeneralHeadBoundary( constant_head, conductance=conductance, concentration=constant_concentration, 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, ) 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="last", save_budget="last") gwf_model["rch"] = imod.mf6.Recharge( rch_rate, concentration=rch_concentration, print_flows=True, save_flows=True ) simulation = imod.mf6.Modflow6Simulation("circle") simulation["flow"] = gwf_model simulation["flow_solver"] = imod.mf6.Solution( modelnames=["flow"], 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, ) .. GENERATED FROM PYTHON SOURCE LINES 191-195 Set the timesteps, we want output each year, so we specify stress periods which last 1 year. However, timesteps of 1 year yield unstable results, so we set ``n_timesteps`` to 10, which sets the amount of timesteps within a stress period. .. GENERATED FROM PYTHON SOURCE LINES 195-200 .. code-block:: Python simtimes = pd.date_range(start="2000-01-01", end="2030-01-01", freq="As") simulation.create_time_discretization(additional_times=simtimes) simulation["time_discretization"]["n_timesteps"] = 10 .. GENERATED FROM PYTHON SOURCE LINES 201-206 Buoyancy -------- Since we are solving a variable density problem, we need to add the buoyancy package. It will use the species "salinity" that we are simulating with a transport model defined below. .. GENERATED FROM PYTHON SOURCE LINES 206-216 .. code-block:: Python slope = (max_density - min_density) / (max_concentration - min_concentration) gwf_model["buoyancy"] = imod.mf6.Buoyancy( reference_density=min_density, modelname=["transport"], reference_concentration=[min_concentration], density_concentration_slope=[slope], species=["salinity"], ) .. GENERATED FROM PYTHON SOURCE LINES 217-225 Add transport model to simulation --------------------------------- The transport model requires the flow field inside the domain computed by the flow model. It also needs the fluxes over the boundary. It uses the same discretization as the flow model. Here we create a transport model for salinity, derive sources and sinks based from the flow model, and tell it to use the same discretization. .. GENERATED FROM PYTHON SOURCE LINES 225-232 .. code-block:: Python transport_model = imod.mf6.GroundwaterTransportModel() transport_model["ssm"] = imod.mf6.SourceSinkMixing.from_flow_model( gwf_model, "salinity" ) transport_model["disv"] = gwf_model["disv"] .. GENERATED FROM PYTHON SOURCE LINES 233-237 Now we define some transport packages for simulating the physical processes of advection, mechanical dispersion, and molecular diffusion dispersion. This example is transient, and the volume available for storage is the porosity, in this case 0.10. .. GENERATED FROM PYTHON SOURCE LINES 237-251 .. code-block:: Python al = 0.001 transport_model["dsp"] = imod.mf6.Dispersion( diffusion_coefficient=1e-4, longitudinal_horizontal=al, transversal_horizontal1=al * 0.1, transversal_vertical=al * 0.01, xt3d_off=False, xt3d_rhs=False, ) transport_model["adv"] = imod.mf6.AdvectionUpstream() transport_model["mst"] = imod.mf6.MobileStorageTransfer(porosity) .. GENERATED FROM PYTHON SOURCE LINES 252-255 Define the maximum concentration as the initial conditions, also output options for the transport model, and assign the transport model to the simulation as well. .. GENERATED FROM PYTHON SOURCE LINES 255-277 .. code-block:: Python transport_model["ic"] = imod.mf6.InitialConditions(start=max_concentration) transport_model["oc"] = imod.mf6.OutputControl( save_concentration="last", save_budget="last" ) simulation["transport"] = transport_model simulation["transport_solver"] = imod.mf6.Solution( modelnames=["transport"], 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, ) .. GENERATED FROM PYTHON SOURCE LINES 278-279 We'll create a new directory in which we will write and run the model. .. GENERATED FROM PYTHON SOURCE LINES 279-283 .. code-block:: Python modeldir = imod.util.temporary_directory() simulation.write(modeldir, binary=False) .. GENERATED FROM PYTHON SOURCE LINES 284-291 Run the model ------------- .. note:: The following lines assume the ``mf6`` executable is available on your PATH. The examples introduction shortly describes how to add it to yours. .. GENERATED FROM PYTHON SOURCE LINES 291-294 .. code-block:: Python simulation.run() .. GENERATED FROM PYTHON SOURCE LINES 295-297 Open the results ---------------- .. GENERATED FROM PYTHON SOURCE LINES 297-301 .. code-block:: Python sim_concentration = simulation.open_concentration().compute() sim_head = simulation.open_head().compute() .. GENERATED FROM PYTHON SOURCE LINES 302-310 Assign coordinates to output ---------------------------- The model output does not feature very useful coordinate values for ``time`` and ``z``, therefore it is best to assign these to the datasets for more understandable plots. First we have to compute values for a z coordinate. The .. GENERATED FROM PYTHON SOURCE LINES 310-315 .. code-block:: Python interfaces = np.concatenate([[top], bottom.values]) z = (interfaces[:-1] + interfaces[1:]) / 2 z .. rst-class:: sphx-glr-script-out .. code-block:: none array([ -5., -15., -25., -35., -45., -55., -65., -75., -85., -95., -105., -115., -125., -135., -145.]) .. GENERATED FROM PYTHON SOURCE LINES 316-317 Assign these new coordinate values to the dataset .. GENERATED FROM PYTHON SOURCE LINES 317-323 .. code-block:: Python coords = {"time": simtimes[1:], "z": ("layer", z)} sim_head = sim_head.assign_coords(**coords) sim_concentration = sim_concentration.assign_coords(**coords) .. GENERATED FROM PYTHON SOURCE LINES 324-329 Visualize the results --------------------- We can quickly and easily visualize the output with the plotting functions provided by xarray and xugrid: .. GENERATED FROM PYTHON SOURCE LINES 329-334 .. code-block:: Python fig, ax = plt.subplots() sim_head.isel(time=-1, layer=0).ugrid.plot(ax=ax) ax.set_aspect(1) .. image-sg:: /examples/mf6/images/sphx_glr_circle_transport_003.png :alt: layer = 1, time = 2030-01-01, z = -5.0 :srcset: /examples/mf6/images/sphx_glr_circle_transport_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 335-337 We can draw a crossection through the center by selecting y=0, for which we can plot the contours as follows: .. GENERATED FROM PYTHON SOURCE LINES 337-342 .. code-block:: Python fig, ax = plt.subplots() sim_concentration.isel(time=-1).ugrid.sel(y=0).plot.contourf( ax=ax, x="mesh2d_x", y="z", cmap="RdYlBu_r" ) .. image-sg:: /examples/mf6/images/sphx_glr_circle_transport_004.png :alt: time = 2030-01-01 :srcset: /examples/mf6/images/sphx_glr_circle_transport_004.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.800 seconds) .. _sphx_glr_download_examples_mf6_circle_transport.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: circle_transport.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: circle_transport.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: circle_transport.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_