ElderΒΆ

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