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
 31# set up constant heads
 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
101# Define starting heads
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)
175m["adv"] = imod.wq.AdvectionTVD(courant=1.0)
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)
181m["pcg"] = imod.wq.PreconditionedConjugateGradientSolver(
182    max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0
183)
184m["gcg"] = imod.wq.GeneralizedConjugateGradientSolver(
185    max_iter=150,
186    inner_iter=30,
187    cclose=1.0e-6,
188    preconditioner="mic",
189    lump_dispersion=True,
190)
191m["oc"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)
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
201# head = imod.idf.open("Hydrocoin/results/head/*.idf")
202# head.plot(yincrease=False)
203# conc = imod.idf.open("Hydrocoin/results/conc/*.idf")
204# conc.plot(levels=range(0, 35, 5), yincrease=False)