Vertical Interface

This 2D examples demonstrates the rotation of an initially vertical interface between fresh and salt water.

For a detailed description of this benchmark, see:

Bakker, M., Oude Essink, G. H. P., & Langevin, C. D. (2004). The rotating movement of three immiscible fluids - A benchmark problem. Journal of Hydrology, 287 (1-4), 270-278. https://doi.org/10.1016/j.jhydrol.2003.10.007

We’ll start with the usual imports

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

import imod

Discretization

We’ll start off by creating a model discretization, since this is a simple conceptual model. The model is a 2D cross-section, hence nrow = 1.

nrow = 1
ncol = 80
nlay = 40

dz = 1.0  # 0.0125
dx = 1.0  # 0.0125
dy = -dx

# Defining tops and bottoms
top1D = xr.DataArray(
    np.arange(nlay * dz, 0.0, -dz), {"layer": np.arange(1, nlay + 1)}, ("layer")
)

bot = top1D - dz
top = nlay * dz

Set up ibound, which sets where active cells are (ibound = 1.0)

bnd = xr.DataArray(
    data=np.full((nlay, nrow, ncol), 1.0),
    coords={
        "y": [0.5],
        "x": np.arange(0.5 * dx, dx * ncol, dx),
        "layer": np.arange(1, 1 + nlay),
        "dx": dx,
        "dy": dy,
    },
    dims=("layer", "y", "x"),
)

fig, ax = plt.subplots()
bnd.plot(y="layer", yincrease=False, ax=ax)

<matplotlib.collections.QuadMesh object at 0x000001A3451BA9D0>

Define the icbund, which sets which cells in the solute transport model are active, inactive or constant.

We just go for active cells everywhere here

icbund = xr.full_like(bnd, 1)

Boundary Conditions

Set the constant heads by specifying a negative value in iboud, that is: bnd[index] = -1`

bnd[31, :, 0] = -1

fig, ax = plt.subplots()
bnd.plot(y="layer", yincrease=False, ax=ax)

<matplotlib.collections.QuadMesh object at 0x000001A34ECDDF50>

Define WEL data

weldata = pd.DataFrame()
weldata["x"] = np.full(1, 0.5 * dx)
weldata["y"] = np.full(1, 0.5)
weldata["q"] = 0.28512  # positive, so it's an injection well

Initial Conditions

Define the starting concentration

sconc = xr.DataArray(
    data=np.full((nlay, nrow, ncol), 0.0),
    coords={
        "y": [0.5],
        "x": np.arange(0.5 * dx, dx * ncol, dx),
        "layer": np.arange(1, nlay + 1),
    },
    dims=("layer", "y", "x"),
)

sconc[:, :, 41:80] = 35.0

Build

Finally, we build the model.

fig, ax = plt.subplots()
sconc.plot(y="layer", yincrease=False, ax=ax)

# Finally, we build the model
m = imod.wq.SeawatModel("VerticalInterface")
m["bas"] = imod.wq.BasicFlow(ibound=bnd, top=top, bottom=bot, starting_head=0.0)
m["lpf"] = imod.wq.LayerPropertyFlow(
    k_horizontal=86.4, k_vertical=86.4, specific_storage=0.0
)
m["btn"] = imod.wq.BasicTransport(
    icbund=icbund, starting_concentration=sconc, porosity=0.1
)
m["adv"] = imod.wq.AdvectionTVD(courant=1.0)
m["dsp"] = imod.wq.Dispersion(longitudinal=0.0, diffusion_coefficient=0.0)
m["vdf"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)
m["wel"] = imod.wq.Well(
    id_name="wel", x=weldata["x"], y=weldata["y"], rate=weldata["q"]
)
m["pcg"] = imod.wq.PreconditionedConjugateGradientSolver(
    max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0
)
m["gcg"] = imod.wq.GeneralizedConjugateGradientSolver(
    max_iter=150,
    inner_iter=30,
    cclose=1.0e-6,
    preconditioner="mic",
    lump_dispersion=True,
)
m["oc"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)
m.create_time_discretization(additional_times=["2000-01-01", "2000-01-02"])

Now we write the model, including runfile:

modeldir = imod.util.temporary_directory()
m.write(modeldir, resultdir_is_workdir=True)

Run

You can run the model using the comand prompt and the iMOD-WQ executable. This is part of the iMOD v5 release, which can be downloaded here: https://oss.deltares.nl/web/imod/download-imod5 . This only works on Windows.

To run your model, open up a command prompt and run the following commands:

cd c:\path\to\modeldir
c:\path\to\imod\folder\iMOD-WQ_V5_3_SVN359_X64R.exe VerticalInterface.run

Note that the version name of your executable might differ.

%% Visualise results —————–

After succesfully running the model, you can plot results as follows:

head = imod.idf.open(modeldir / "results/head/*.idf")

fig, ax = plt.subplots()
head.plot(yincrease=False, ax=ax)

conc = imod.idf.open(modeldir / "results/conc/*.idf")

fig, ax = plt.subplots()
conc.plot(levels=range(0, 35, 5), yincrease=False, ax=ax)