Freshwater lensΒΆ

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