TWRI Modflow 6 example 1ΒΆ

  1import subprocess
  2from pathlib import Path
  3
  4import numpy as np
  5import xarray as xr
  6
  7import imod
  8
  9nlay = 3
 10nrow = 15
 11ncol = 15
 12shape = (nlay, nrow, ncol)
 13
 14dx = 5000.0
 15dy = -5000.0
 16xmin = 0.0
 17xmax = dx * ncol
 18ymin = 0.0
 19ymax = abs(dy) * nrow
 20dims = ("layer", "y", "x")
 21
 22layer = np.array([1, 2, 3])
 23y = np.arange(ymax, ymin, dy) + 0.5 * dy
 24x = np.arange(xmin, xmax, dx) + 0.5 * dx
 25coords = {"layer": layer, "y": y, "x": x}
 26
 27# Discretization data
 28idomain = xr.DataArray(np.ones(shape), coords=coords, dims=dims)
 29bottom = xr.DataArray([-200.0, -300.0, -450.0], {"layer": layer}, ("layer",))
 30
 31# Constant head
 32head = xr.full_like(idomain, np.nan).sel(layer=[1, 2])
 33head[...] = np.nan
 34head[..., 0] = 0.0
 35
 36# Drainage
 37elevation = xr.full_like(idomain.sel(layer=1), np.nan)
 38conductance = xr.full_like(idomain.sel(layer=1), np.nan)
 39elevation[7, 1:10] = np.array([0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0])
 40conductance[7, 1:10] = 1.0
 41
 42# Node properties
 43icelltype = xr.DataArray([1, 0, 0], {"layer": layer}, ("layer",))
 44k = xr.DataArray([1.0e-3, 1.0e-4, 2.0e-4], {"layer": layer}, ("layer",))
 45k33 = xr.DataArray([2.0e-8, 2.0e-8, 2.0e-8], {"layer": layer}, ("layer",))
 46
 47# Recharge
 48rch_rate = xr.full_like(idomain.sel(layer=1), 3.0e-8)
 49
 50# Well
 51layer = [3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
 52row = [5, 4, 6, 9, 9, 9, 9, 11, 11, 11, 11, 13, 13, 13, 13]
 53column = [11, 6, 12, 8, 10, 12, 14, 8, 10, 12, 14, 8, 10, 12, 14]
 54rate = [-5.0] * 15
 55
 56# Create and fill the groundwater model.
 57gwf_model = imod.mf6.GroundwaterFlowModel()
 58gwf_model["dis"] = imod.mf6.StructuredDiscretization(
 59    top=200.0, bottom=bottom, idomain=idomain
 60)
 61gwf_model["chd"] = imod.mf6.ConstantHead(
 62    head, print_input=True, print_flows=True, save_flows=True
 63)
 64gwf_model["drn"] = imod.mf6.Drainage(
 65    elevation=elevation,
 66    conductance=conductance,
 67    print_input=True,
 68    print_flows=True,
 69    save_flows=True,
 70)
 71gwf_model["ic"] = imod.mf6.InitialConditions(head=0.0)
 72gwf_model["npf"] = imod.mf6.NodePropertyFlow(
 73    icelltype=icelltype,
 74    k=k,
 75    k33=k33,
 76    variable_vertical_conductance=True,
 77    dewatered=True,
 78    perched=True,
 79    save_flows=True,
 80)
 81gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all")
 82gwf_model["rch"] = imod.mf6.Recharge(rch_rate)
 83gwf_model["wel"] = imod.mf6.Well(
 84    layer=layer,
 85    row=row,
 86    column=column,
 87    rate=rate,
 88    print_input=True,
 89    print_flows=True,
 90    save_flows=True,
 91)
 92
 93# Attach it to a simulation
 94simulation = imod.mf6.Modflow6Simulation("ex01-twri")
 95simulation["GWF_1"] = gwf_model
 96# Define solver settings
 97simulation["solver"] = imod.mf6.Solution(
 98    print_option="summary",
 99    csv_output=False,
100    no_ptc=True,
101    outer_hclose=1.0e-4,
102    outer_maximum=500,
103    under_relaxation=None,
104    inner_hclose=1.0e-4,
105    inner_rclose=0.001,
106    inner_maximum=100,
107    linear_acceleration="cg",
108    scaling_method=None,
109    reordering_method=None,
110    relaxation_factor=0.97,
111)
112# Collect time discretization
113simulation.time_discretization(times=["2000-01-01", "2000-01-02"])
114
115modeldir = Path("ex01-twri")
116simulation.write(modeldir)
117with imod.util.cd(modeldir):
118    p = subprocess.run("mf6", check=True, capture_output=True, text=True)
119    assert p.stdout.endswith("Normal termination of simulation.\n")
120
121head = imod.mf6.open_hds(modeldir / "GWF_1/GWF_1.hds", modeldir / "GWF_1/dis.dis.grb")