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
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)
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)
83    max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0
84)
86    max_iter=150,
87    inner_iter=30,
88    cclose=1.0e-6,
89    preconditioner="mic",
90    lump_dispersion=True,
91)