Vertical interfaceΒΆ

  1import numpy as np
  2import pandas as pd
  3import xarray as xr
  4
  5import imod
  6
  7# Discretization
  8nrow = 1
  9ncol = 80
 10nlay = 40
 11
 12dz = 1.0  # 0.0125
 13dx = 1.0  # 0.0125
 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# Defining constant heads
 32bnd[31, :, 0] = -1
 33bnd.plot(y="layer", yincrease=False)
 34
 35# Defining tops and bottoms
 36top1D = xr.DataArray(
 37    np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, nlay + 1)}, ("layer")
 38)
 39
 40bot = top1D - dz
 41top = nlay * dz
 42
 43# Define WEL data
 44weldata = pd.DataFrame()
 45weldata["x"] = np.full(1, 0.5 * dx)
 46weldata["y"] = np.full(1, 0.5)
 47weldata["q"] = 0.28512  # positive, so it's an injection well
 48
 49# Define icbund
 50icbund = xr.full_like(bnd, 1)
 51
 52# Define starting concentrations
 53sconc = xr.DataArray(
 54    data=np.full((nlay, nrow, ncol), 0.0),
 55    coords={
 56        "y": [0.5],
 57        "x": np.arange(0.5 * dx, dx * ncol, dx),
 58        "layer": np.arange(1, nlay + 1),
 59    },
 60    dims=("layer", "y", "x"),
 61)
 62
 63sconc[:, :, 41:80] = 35.0
 64sconc.plot(y="layer", yincrease=False)
 65
 66# Finally, we build the model
 67m = imod.wq.SeawatModel("VerticalInterface")
 68m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=top, bottom=bot, starting_head=0.0)
 69m["lpf"] = imod.wq.LayerPropertyFlow(
 70    k_horizontal=86.4, k_vertical=86.4, specific_storage=0.0
 71)
 72m["btn"] = imod.wq.BasicTransport(
 73    icbund=icbund, starting_concentration=sconc, porosity=0.1
 74)
 75m["adv"] = imod.wq.AdvectionTVD(courant=1.0)
 76m["dsp"] = imod.wq.Dispersion(longitudinal=0.0, diffusion_coefficient=0.0)
 77m["vdf"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)
 78m["wel"] = imod.wq.Well(
 79    id_name="wel", x=weldata["x"], y=weldata["y"], rate=weldata["q"]
 80)
 81m["pcg"] = imod.wq.PreconditionedConjugateGradientSolver(
 82    max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0
 83)
 84m["gcg"] = imod.wq.GeneralizedConjugateGradientSolver(
 85    max_iter=150,
 86    inner_iter=30,
 87    cclose=1.0e-6,
 88    preconditioner="mic",
 89    lump_dispersion=True,
 90)
 91m["oc"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)
 92m.time_discretization(times=["2000-01-01", "2000-01-02"])
 93
 94# Now we write the model, including runfile:
 95m.write()
 96
 97
 98# You can run the model using the comand prompt and the iMOD SEAWAT executable
 99
100# Results
101
102# head = imod.idf.open("VerticalInterface/results/head/*.idf")
103# head.plot(yincrease=False)
104# conc = imod.idf.open("VerticalInterface/results/conc/*.idf")
105# conc.plot(levels=range(0, 35, 5), yincrease=False)