Note
Go to the end to download the full example code.
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).
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