# Hydro coinΒΆ

```  1import numpy as np
2import pandas as pd
3import xarray as xr
4
5import imod
6
7# Discretization
8nrow = 1  # number of rows
9ncol = 45  # number of columns
10nlay = 76  # number of layers
11
12dz = 4.0
13dx = 20.0
14dy = -dx
15
16# setup ibound
17bnd = xr.DataArray(
18    data=np.full((nlay, nrow, ncol), 1.0),
19    coords={
20        "y": [0.5],
21        "x": np.arange(0.5 * dx, dx * ncol, dx),
22        "layer": np.arange(1, 1 + nlay),
23        "dx": dx,
24        "dy": dy,
25    },
26    dims=("layer", "y", "x"),
27)
28
29bnd.plot()
30
32bnd[0, :, :] = -1
33bnd[75, :, 0:15] = 0.0
34bnd[75, :, 30:45] = 0.0
35bnd.plot(y="layer", yincrease=False)
36
37# setup tops and bottoms
38top1D = xr.DataArray(
39    np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, nlay + 1)}, ("layer")
40)
41
42bot = top1D - dz
43
44# Define WEL data, need to define the x, y, and pumping rate (q)
45weldata = pd.DataFrame()
46weldata["x"] = np.full(1, 0.5 * dx)
47weldata["y"] = np.full(1, 0.5)
48weldata["q"] = 0.28512  # positive, so it's an injection well
49
50
51# Define the starting concentrations
52sconc = xr.DataArray(
53    data=np.full((nlay, nrow, ncol), 0.0),
54    coords={
55        "y": [0.5],
56        "x": np.arange(0.5 * dx, dx * ncol, dx),
57        "layer": np.arange(1, nlay + 1),
58        "dx": dx,
59        "dy": dy,
60    },
61    dims=("layer", "y", "x"),
62)
63
64sconc[75, :, 15:30] = 280.0
65sconc.plot(y="layer", yincrease=False)
66
67# Define the icbund
68icbund = xr.DataArray(
69    data=np.full((nlay, nrow, ncol), 1.0),
70    coords={
71        "y": [0.5],
72        "x": np.arange(0.5 * dx, dx * ncol, dx),
73        "layer": np.arange(1, nlay + 1),
74        "dx": dx,
75        "dy": dy,
76    },
77    dims=("layer", "y", "x"),
78)
79
80icbund[75, :, 0:15] = 0.0
81icbund[75, :, 30:45] = 0.0
82icbund[75, :, 15:30] = -1.0
83icbund.plot(y="layer", yincrease=False)
84
85# Define horizontal hydraulic conductivity
86khv = xr.DataArray(
87    data=np.full((nlay, nrow, ncol), 0.847584),
88    coords={
89        "y": [0.5],
90        "x": np.arange(0.5 * dx, dx * ncol, dx),
91        "layer": np.arange(1, nlay + 1),
92        "dx": dx,
93        "dy": dy,
94    },
95    dims=("layer", "y", "x"),
96)
97
98khv[75, :, 15:30] = 0.0008475
99khv.plot(y="layer", yincrease=False)
100
102shd = xr.DataArray(
103    data=np.full((nlay, nrow, ncol), 0.0),
104    coords={
105        "y": [0.5],
106        "x": np.arange(0.5 * dx, dx * ncol, dx),
107        "layer": np.arange(1, nlay + 1),
108        "dx": dx,
109        "dy": dy,
110    },
111    dims=("layer", "y", "x"),
112)
113
114shd[0, :, :] = np.array(
115    [
116        10,
117        9.772727273,
118        9.545454545,
119        9.318181818,
120        9.090909091,
121        8.863636364,
122        8.636363636,
123        8.409090909,
124        8.181818182,
125        7.954545455,
126        7.727272727,
127        7.5,
128        7.272727273,
129        7.045454545,
130        6.818181818,
131        6.590909091,
132        6.363636364,
133        6.136363636,
134        5.909090909,
135        5.681818182,
136        5.454545455,
137        5.227272727,
138        5,
139        4.772727273,
140        4.545454545,
141        4.318181818,
142        4.090909091,
143        3.863636364,
144        3.636363636,
145        3.409090909,
146        3.181818182,
147        2.954545455,
148        2.727272727,
149        2.5,
150        2.272727273,
151        2.045454545,
152        1.818181818,
153        1.590909091,
154        1.363636364,
155        1.136363636,
156        0.909090909,
157        0.681818182,
158        0.454545455,
159        0.227272727,
160        0.00,
161    ]
162)
163shd.plot(y="layer", yincrease=False)
164
165
166# Finally, we build the model.
167m = imod.wq.SeawatModel("Hydrocoin")
168m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=304.0, bottom=bot, starting_head=shd)
169m["lpf"] = imod.wq.LayerPropertyFlow(
170    k_horizontal=khv, k_vertical=khv, specific_storage=0.0
171)
172m["btn"] = imod.wq.BasicTransport(
173    icbund=icbund, starting_concentration=sconc, porosity=0.2
174)
176m["dsp"] = imod.wq.Dispersion(longitudinal=20.0, diffusion_coefficient=0.0)
177m["vdf"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)
178m["wel"] = imod.wq.Well(
179    id_name="wel", x=weldata["x"], y=weldata["y"], rate=weldata["q"]
180)
182    max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0
183)
185    max_iter=150,
186    inner_iter=30,
187    cclose=1.0e-6,
188    preconditioner="mic",
189    lump_dispersion=True,
190)
192m.time_discretization(times=["2000-01-01T00:00", "2010-01-01T00:00"])
193
194
195# Now we write the model, including runfile:
196m.write()
197
198# You can run the model using the comand prompt and the iMOD SEAWAT executable
199
200## Visualise results