Salt water pocketΒΆ

 1import numpy as np
 2import pandas as pd
 3import xarray as xr
 4
 5import imod
 6
 7# Discretization
 8nrow = 1  # number of rows
 9ncol = 80  # number of column
10nlay = 40  # number of layers
11
12dz = 1.0  # 0.0125
13dx = 1.0  # 0.0125
14dy = -dx
15
16
17# setup ibound
18bnd = xr.DataArray(
19    data=np.full((nlay, nrow, ncol), 1.0),
20    coords={
21        "y": [0.5],
22        "x": np.arange(0.5 * dx, dx * ncol, dx),
23        "layer": np.arange(1, 1 + nlay),
24        "dx": dx,
25        "dy": dy,
26    },
27    dims=("layer", "y", "x"),
28)
29
30bnd.plot()
31
32# set constant heads
33bnd[21, :, 0] = -1
34
35# setup tops and bottoms
36
37top1D = xr.DataArray(
38    np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, nlay + 1)}, ("layer")
39)
40
41bot = top1D - dz
42
43# Defining the starting concentrations
44sconc = xr.DataArray(
45    data=np.full((nlay, nrow, ncol), 0.0),
46    coords={
47        "y": [0.5],
48        "x": np.arange(0.5 * dx, dx * ncol, dx),
49        "layer": np.arange(1, nlay + 1),
50    },
51    dims=("layer", "y", "x"),
52)
53
54sconc[16:24, :, 41:80] = 35.0
55sconc.plot(y="layer", yincrease=False)
56
57# Finally, we build the model.
58
59m = imod.wq.SeawatModel("SaltwaterPocket")
60m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=40, bottom=bot, starting_head=0.0)
61m["lpf"] = imod.wq.LayerPropertyFlow(
62    k_horizontal=86.4, k_vertical=86.4, specific_storage=0.0
63)
64m["btn"] = imod.wq.BasicTransport(
65    icbund=bnd, starting_concentration=sconc, porosity=0.1
66)
67m["adv"] = imod.wq.AdvectionTVD(courant=1.0)
68m["dsp"] = imod.wq.Dispersion(longitudinal=0.001, diffusion_coefficient=0.0000864)
69m["vdf"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)
70m["wel"] = imod.wq.Well(id_name="wel", x=0.5 * dx, y=0.5, rate=0.28512)
71m["pcg"] = imod.wq.PreconditionedConjugateGradientSolver(
72    max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0
73)
74m["gcg"] = imod.wq.GeneralizedConjugateGradientSolver(
75    max_iter=150,
76    inner_iter=30,
77    cclose=1.0e-6,
78    preconditioner="mic",
79    lump_dispersion=True,
80)
81m["oc"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)
82m.time_discretization(times=["2000-01-01T00:00", "2000-01-05T01:00"])
83
84
85# Now we write the model, including runfile:
86
87m.write()
88
89
90# You can run the model using the comand prompt and the iMOD SEAWAT executable
91
92# Results
93
94# head = imod.idf.open("SaltwaterPocket/results/head/*.idf")
95# head.plot(yincrease=False)
96# conc = imod.idf.open("SaltwaterPocket/results/conc/*.idf")
97# conc.plot(levels=range(0, 35, 5), yincrease=False)