.. 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 <sphx_glr_download_examples_mf6_Henry.py>`
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 MODFLOW 6 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 do 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 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 <mf6-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
<matplotlib.contour.QuadContourSet object at 0x0000016F3FDA2F60>
.. 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
<matplotlib.contour.QuadContourSet object at 0x0000016F3FD4D7C0>
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 58.243 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 <Henry.ipynb>`
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: Henry.py <Henry.py>`
.. container:: sphx-glr-download sphx-glr-download-zip
:download:`Download zipped: Henry.zip <Henry.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_