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