1import numpy as np
2import pandas as pd
3import xarray as xr
4
5import imod
6
7# Discretization
8nrow = 1
9ncol = 100
10nlay = 50
11
12dz = 1.0
13dx = 1.0
14dy = -dx
15
16# scale parameters with discretization
17qscaled = 0.03 * (dz * abs(dy))
18
19# Fresh water injection with well
20# Add the arguments as a list, so pandas doesn't complain about having to set
21# an index.
22weldata = pd.DataFrame()
23weldata["x"] = [0.5]
24weldata["y"] = [0.5]
25weldata["q"] = [qscaled]
26
27# Setup ibound
28bnd = xr.DataArray(
29 data=np.full((nlay, nrow, ncol), 1.0),
30 coords={
31 "y": [0.5],
32 "x": np.arange(0.5 * dx, dx * ncol, dx),
33 "layer": np.arange(1, 1 + nlay),
34 "dx": dx,
35 "dy": dy,
36 },
37 dims=("layer", "y", "x"),
38)
39
40top1D = xr.DataArray(
41 np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, 1 + nlay)}, ("layer")
42)
43bot = top1D - 1.0
44# We define constant head here, after generating the tops, or we'd end up with negative top values
45bnd[:, :, -1] = -1
46
47# Fill model
48m = imod.wq.SeawatModel("HenryCase")
49m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=50.0, bottom=bot, starting_head=1.0)
50m["lpf"] = imod.wq.LayerPropertyFlow(
51 k_horizontal=10.0, k_vertical=10.0, specific_storage=0.0
52)
53m["btn"] = imod.wq.BasicTransport(
54 icbund=bnd, starting_concentration=35.0, porosity=0.35
55)
56m["adv"] = imod.wq.AdvectionTVD(courant=1.0)
57m["dsp"] = imod.wq.Dispersion(longitudinal=0.1, diffusion_coefficient=1.0e-9)
58m["vdf"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)
59m["wel"] = imod.wq.Well(
60 id_name="well", x=weldata["x"], y=weldata["y"], rate=weldata["q"]
61)
62m["pcg"] = imod.wq.PreconditionedConjugateGradientSolver(
63 max_iter=150, inner_iter=30, hclose=0.0001, rclose=1.0, relax=0.98, damp=1.0
64)
65m["gcg"] = imod.wq.GeneralizedConjugateGradientSolver(
66 max_iter=150,
67 inner_iter=30,
68 cclose=1.0e-6,
69 preconditioner="mic",
70 lump_dispersion=True,
71)
72m["oc"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)
73m.time_discretization(times=pd.date_range("2000-01-01", "2001-01-01", freq="M"))
74m.write("HenryCase")