.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\mf6\Henry.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_Henry.py: Henry ===== This example illustrates how to setup a variable density groundwater flow and transport model using the ``imod`` package and associated packages. In overview, we'll set the following steps: * Create a suitable 2d (x, z) grid. * Create a groundwater flow model, with variable density. * Create a solute transport model. * Combine these models into a single MODFLOW6 simulation. * Write to modflow6 files. * Run the model. * Open the results back into xarray DataArrays. * Visualize the results. We are simulating the Henry problem, although not the original one but the one outlined in the MODFLOW 6 manual (jupyter notebooks example 51). This is the modified Henry problem with half the inflow rate. The domain is a vertically oriented two dimensional rectangle, which is 2 m long and 1 m high. Water flows in over the left boundary with a fixed rate, which is represented by a Well package. The right boundary is in direct contact with hydrostatic seawater with a density of 1025 kg m:sup:`-3`. This is represented by a General Head Boundary package. .. GENERATED FROM PYTHON SOURCE LINES 29-31 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 33-35 We'll start with the usual imports. As this is a simple (synthetic) structured model, we can make due with few packages. .. GENERATED FROM PYTHON SOURCE LINES 35-42 .. code-block:: Python import numpy as np import pandas as pd import xarray as xr import imod .. GENERATED FROM PYTHON SOURCE LINES 43-45 We'll start by defining the (vertical) rectangular domain and the physical parameters of the model. .. GENERATED FROM PYTHON SOURCE LINES 45-90 .. code-block:: Python nlay = 40 nrow = 1 ncol = 80 shape = (nlay, nrow, ncol) total_flux = 5.7024 # m3/d k = 864.0 # m/d porosity = 0.35 max_concentration = 35.0 min_concentration = 0.0 max_density = 1025.0 min_density = 1000.0 diffusion_coefficient = 0.57024 longitudinal_horizontal = 0.1 transversal_horizontal1 = 0.01 # Time start_date = pd.to_datetime("2020-01-01") duration = pd.to_timedelta("0.5d") # Domain size xmax = 2.0 xmin = 0.0 dx = (xmax - xmin) / ncol zmin = 0.0 zmax = 1.0 dz = (zmax - zmin) / nlay x = np.arange(xmin, xmax, dx) + 0.5 * dx y = np.array([0.5]) layer = np.arange(1, 41, 1) dy = -1.0 coords = {"layer": layer, "y": y, "x": x, "dy": dy, "dx": dx} dims = ("layer", "y", "x") idomain = xr.DataArray(np.ones(shape, dtype=int), coords=coords, dims=dims) top = 1.0 bottom = xr.DataArray( np.arange(zmin, zmax, dz)[::-1], {"layer": layer}, ("layer",), ) .. GENERATED FROM PYTHON SOURCE LINES 91-93 Now make the flow model. We'll start with the non-boundary condition packages. .. GENERATED FROM PYTHON SOURCE LINES 93-111 .. code-block:: Python gwf_model = imod.mf6.GroundwaterFlowModel() gwf_model["dis"] = imod.mf6.StructuredDiscretization( top=top, bottom=bottom, idomain=idomain ) gwf_model["npf"] = imod.mf6.NodePropertyFlow( icelltype=0, k=k, ) gwf_model["sto"] = imod.mf6.SpecificStorage( specific_storage=1.0e-4, specific_yield=0.15, transient=False, convertible=0, ) gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0) gwf_model["oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last") .. GENERATED FROM PYTHON SOURCE LINES 112-113 Now let's make the constant head boundary condition. .. GENERATED FROM PYTHON SOURCE LINES 113-134 .. code-block:: Python ghb_head = xr.ones_like(idomain, dtype=float) ghb_head[:, :, :-1] = np.nan ghb_conc = xr.full_like(idomain, max_concentration, dtype=float) ghb_conc[:, :, :-1] = np.nan ghb_conc = ghb_conc.expand_dims(species=["salinity"]) conductance = xr.full_like(idomain, 864.0 * 2.0, dtype=float) conductance[:, :, :-1] = np.nan gwf_model["right_boundary"] = imod.mf6.GeneralHeadBoundary( head=ghb_head, conductance=conductance, concentration=ghb_conc, concentration_boundary_type="AUX", print_input=True, print_flows=True, save_flows=True, ) .. GENERATED FROM PYTHON SOURCE LINES 135-136 ... and the constant flux condition. .. GENERATED FROM PYTHON SOURCE LINES 136-157 .. code-block:: Python from imod.prepare.layer import create_layered_top screen_top = create_layered_top(bottom, top) flux_concentration = xr.DataArray( data=np.full((1, nlay), min_concentration), dims=["species", "index"], coords={"species": ["salinity"], "index": layer}, ) gwf_model["left_boundary"] = imod.mf6.Well( x=np.full_like(layer, xmin, dtype=float), y=np.full_like(layer, y[0], dtype=float), screen_top=screen_top.values, screen_bottom=bottom.values, rate=np.full_like(layer, 0.5 * (total_flux / nlay), dtype=float), concentration=flux_concentration, concentration_boundary_type="AUX", minimum_thickness=0.002, ) .. GENERATED FROM PYTHON SOURCE LINES 158-161 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 161-171 .. 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 172-176 Now let's make the transport model. It contains the standard packages of storage, dispersion and advection as well as initial condiations and output control. Sinks and sources are automatically determined based on packages provided in the flow model. .. GENERATED FROM PYTHON SOURCE LINES 176-195 .. code-block:: Python gwt_model = imod.mf6.GroundwaterTransportModel() gwt_model["ssm"] = imod.mf6.SourceSinkMixing.from_flow_model(gwf_model, "salinity") gwt_model["adv"] = imod.mf6.AdvectionTVD() gwt_model["dsp"] = imod.mf6.Dispersion( diffusion_coefficient=0.57024, longitudinal_horizontal=0.1, transversal_horizontal1=0.01, xt3d_off=False, xt3d_rhs=False, ) gwt_model["mst"] = imod.mf6.MobileStorageTransfer( porosity=porosity, ) gwt_model["ic"] = imod.mf6.InitialConditions(start=max_concentration) gwt_model["oc"] = imod.mf6.OutputControl(save_concentration="last", save_budget="last") gwt_model["dis"] = gwf_model["dis"] .. GENERATED FROM PYTHON SOURCE LINES 196-197 now let's define a simulation using the flow and transport models. .. GENERATED FROM PYTHON SOURCE LINES 197-204 .. code-block:: Python # Attach it to a simulation simulation = imod.mf6.Modflow6Simulation("henry") simulation["flow"] = gwf_model simulation["transport"] = gwt_model .. GENERATED FROM PYTHON SOURCE LINES 205-209 Define solver settings. We need to define separate solutions for the flow and transport models. In this case, we'll use the same settings, but generally convergence settings should differ: the transport model has very different units from the flow model. .. GENERATED FROM PYTHON SOURCE LINES 209-242 .. code-block:: Python simulation["flow_solver"] = imod.mf6.Solution( modelnames=["flow"], print_option="summary", outer_dvclose=1.0e-6, outer_maximum=500, under_relaxation=None, inner_dvclose=1.0e-6, inner_rclose=1.0e-10, inner_maximum=100, linear_acceleration="bicgstab", scaling_method=None, reordering_method=None, relaxation_factor=0.97, ) simulation["transport_solver"] = imod.mf6.Solution( modelnames=["transport"], print_option="summary", outer_dvclose=1.0e-6, outer_maximum=500, under_relaxation=None, inner_dvclose=1.0e-6, inner_rclose=1.0e-10, inner_maximum=100, linear_acceleration="bicgstab", scaling_method=None, reordering_method=None, relaxation_factor=0.97, ) # Collect time discretization times = [start_date, start_date + duration] simulation.create_time_discretization(additional_times=times) .. GENERATED FROM PYTHON SOURCE LINES 243-244 Increase the number of time steps for the single stress period: .. GENERATED FROM PYTHON SOURCE LINES 244-246 .. code-block:: Python simulation["time_discretization"]["n_timesteps"] = 500 .. GENERATED FROM PYTHON SOURCE LINES 247-248 We'll create a new directory in which we will write and run the model. .. GENERATED FROM PYTHON SOURCE LINES 248-251 .. code-block:: Python modeldir = imod.util.temporary_directory() simulation.write(modeldir, binary=False) .. GENERATED FROM PYTHON SOURCE LINES 252-262 Run the model ------------- This takes about 20 seconds. .. 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 262-265 .. code-block:: Python simulation.run() .. GENERATED FROM PYTHON SOURCE LINES 266-270 Open the results ---------------- We'll open the head and concentration files. .. GENERATED FROM PYTHON SOURCE LINES 270-275 .. code-block:: Python head = simulation.open_head() conc = simulation.open_concentration() .. GENERATED FROM PYTHON SOURCE LINES 276-278 Visualize the results --------------------- .. GENERATED FROM PYTHON SOURCE LINES 278-281 .. code-block:: Python head.isel(y=0, time=-1).plot.contourf(yincrease=False) .. image-sg:: /examples/mf6/images/sphx_glr_Henry_001.png :alt: y = 0.5, dx = 0.025, dy = -1.0, time = 0.5 :srcset: /examples/mf6/images/sphx_glr_Henry_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 282-284 We can check the concentration to see that a fresh-saline interface has been formed: .. GENERATED FROM PYTHON SOURCE LINES 284-287 .. code-block:: Python conc.isel(y=0, time=-1).plot.contourf(yincrease=False, cmap="RdYlBu_r") .. image-sg:: /examples/mf6/images/sphx_glr_Henry_002.png :alt: y = 0.5, dx = 0.025, dy = -1.0, time = 0.5 :srcset: /examples/mf6/images/sphx_glr_Henry_002.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.098 seconds) .. _sphx_glr_download_examples_mf6_Henry.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Henry.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Henry.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Henry.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_